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Periodic Motions of a Rectangular Wing 
Moving at Supersonic Speed 


H. J. STEWART* ano TING-YI LIt 
California Institute of Technology 


(1) SumMMARY 


The determination of the air forces and moments on an oscillat- 
ing airfoil in a main stream of uniform supersonic speed is im- 
portant because of its vital role in the prediction of the flutter 
and aerodynamic instability characteristics of high-speed air- 
craft and supersonic missiles. In the past, calculations concern- 
ing this problem were performed mostly on the basis of an analysis 
of a two-dimensional wing of infinite span. A satisfactory theory 
of the corresponding problem for a three-dimensional wing with 
finite span is not yet available. In reference 1, the present au 
thors developed some general integral expressions for the dis 
turbance velocity potential of a three-dimensional oscillating 
wing with both purely and mixed supersonic boundary condi 
tions. The method used was essentially that of Evvard.? 
These general results are applied in the present paper to deter 
mine the lift and the moment (due to the lift) acting upon a thin 
oscillating wing of rectangular plan form moving at a supersonic 
speed. Detailed analysis in the case of vertical oscillations is 
presented. Progress in the analysis of pitching oscillations is also 
indicated. 

A fundamentally different method of approach to the oscillating 
rectangular wing problem has been presented by Miles.5 The 


numerical results presented here are consistent with results (not 
yet published) obtained by his method. A method similar to 
that used here is being used by C. C. Chang to compute an ap- 
proximate solution of rectangular wing problems. 


SYMBOLS AND NOTATIONS 


Cartesian coordinates 

running Cartesian coordinates 

time variable 

disturbance velocity potential or ele- 
mentary oscillating source potential 

velocity potential on a point on the wing 
surface 

speed of sound in the free stream 

free-stream speed in the positive x-direc- 
tion 

free-stream Mach Number 


Presented at the Aerodynamics Session, Eighteenth Annual 
Meeting, I.A.S., New York, January 23-26, 1950. 


w(t, n) exp (ivt) 


M? — 1 = §? 

frequency of the periodic motions of the 
wing 

upwash velocity at a point (£, 7) on the 
top surface of the wing at an instant ¢ 

characteristic coordinates 

rectangular wing chord 

rectangular wing span 

rectangular wing aspect ratio 

rectangular wing angle of attack 


reduced frequency 


disturbance pressure 

free-stream density 

lift, moment coefficients in the region II 
of the rectangular wing (Fig. 3) 

lift, moment coefficients in the region I 
of the rectangular wing (Fig. 3) 

lift, moment coefficients in the region 
II’ of the rectangular wing (Fig. 3) 

lift coefficient of the wing, vertical oscilla- 
tion 

moment coefficient of the wing, vertical 
oscillation 

lift coefficient of the wing, pitching os 

» cillation 

= moment coefficient of the wing, pitching 


oscillation 


(III GENERAL INTEGRAL SOLUTIONS! 


Pie A THIN OSCILLATING WING at a small angle of 


attack in a main stream of uniform supersonic speed 


along the x-axis, the basic linearized equation of the 
disturbance velocity potential is 


1 21 
. Pu + . Pri + (\? —_ liz _ Pyy + Pez (1) 
a- a~ 


* Professor of Aeronautics. An elementary solution that has a simple harmonic time 


t Graduate Student. dependence is 
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Fic. 1. Arbitrary lifting wing in a supersonic flow. 


exp {iv[t — (Ux/B?a?)] | 
[x? — B(y? + z?)|'/* 
.? ... a mit «a 
cos , —— [x? — B7(y? + 27)] { (2) 
Ba 
This solution, in the limit vy — 0, reduces to the well- 
known steady-state supersonic source potential. For 
v ~ 0, Eq. (2) may be considered as defining a super- 


G(x, y, 2, t) = 


sonic oscillating source. 

The representation of the steady, two- or three- 
dimensional flow past a wing by distributions of steady 
supersonic sources has been discussed by many authors. 
Probably the most interesting general treatment of this 
problem was given by Evvard.** Source superposi- 
tion methods have also been used for representing super- 
sonic nonstationary motions by Garrick and Rubinow* 
and by Miles.’ Recently, Evvard* has studied the 
nonstationary problems by the same techniques that 
he had previously applied to steady-flow problems. 
Evvard’s approach has been applied by the present 
authors! specifically to oscillating problems, and simple 
results, closely paralleling Evvard’s steady-state theory, 
were obtained. The purpose of the present work is to 
apply this general theory to a particular example, the 
rectangular flat plate wing. 

The results of the general theory of reference | can 
readily be applied to the computation of the perturba- 
tion velocity potential on the surface of an arbitrary 
flat plate wing such as that shown in Fig. |. The lead- 
ing edge in Fig. 1 is divided into three segments. The 
segment OO’ is a supersonic leading edge and the two 
segments, OA and O’B, are subsonic leading edges. It 
was assumed that the entire trailing edge AB is super- 
sonic and that the two subsonic leading edges are inde- 
pendent (i.e., that any Mach wave emitted from a point 
on OA will not intersect O’B and vice versa). The 
present analysis deals exclusively with the lift problem. 
For this case the wing is assumed to lie near the plane 
z = 0, and the perturbation velocity potential is an odd 
function of z. 

It is shown in reference | that for an oscillating wing 
it is possible to replace the wing by a distribution of 
oscillating sources in the plane z = 0; therefore, the 


complete perturbation potential is 


P(x, y, 2, ¢) = [eu A(é, n) Ox — & vy — 7, 3, t) dtdy 
Ss gn 


where the elementary potential ¢ is given by Eq. (2) ang 

A(&, n) is the source strength per unit area. In gen. 

eral, the integration area S is that portion of the plane 

z = 0 intercepted by the forward Mach cone from the 

point (x, y, 2). In this case, the boundary condition 

satisfied by Eq. (3) on the top side of the plane ¢ = ()jg 
lim (08/02) = —mA(x, y) exp (i) { 
z= +0 


If the boundary condition on the wing requires that 


lim (O&/0z) = w(x, y) exp (a1) 5) 
z= +0 
i.e., if the wing is oscillating in some manner—the te 
quired source strength is, thus, 


A(x, y) = — (1/r)w(s, y) (6) 


For a point on the top surface of the wing in the 
purely supersonic region, such as P in Fig. 1, Eqs. (3 
and (6) represent the complete solution. Since the 
disturbance in the portion of the Mach cone from P 
ahead of the wing is zero, the region of integration S is 
bounded by the two characteristics PP; and PP, and 
the leading edge P,P». 

For a surface point such as Q in the mixed supersonic 
region influenced by the subsonic leading edge OA, the 
analysis is more complex, since there is an upwash of 
unknown magnitude in the region bounded by the 
leading edge OA and the Mach wave Ov. In order 
to apply the formula for the perturbation potential 
at point P to point Q, it would be necessary to 
determine the upwash velocity in the region bounded 
by OA and Ov and thus to determine [by Eq. (6)] the 
appropriate source distribution for this region. For- 
tunately, it is not necessary to compute this upwash 
velocity explicitly. It is shown in reference | that the 
contribution of this upwash field (i.e., the effect of the 
bottom of the wing) to the perturbation potential at 
point Q on the top surface just cancels the contribution 
to the perturbation potential at Q due to the sources in 
the region of the wing forward of the Mach line Q:Q:. 
The velocity potential at point Q is thus given by Eqs. 
(3) and (6), where the region of integration S is bounded 
by the characteristics QQ;, QQ:, and Q,Q» and the 
leading edge Q.Q3. 

For a point influenced by only the left-hand sub- 
sonic leading edge, the integration area is similarly 
determined. A point in the triangle DEF is influenced 
by both wing tips and portions of the wing surface on 
both right- and left-hand sides must be omitted in de- 
termining the appropriate area S. 

The appropriate integration area S for an oscillating 
wing is thus identical to that determined by Evvard? for 
stationary wings. It is shown in reference | that, 
while this result is true for steady-flow and for oscil- 
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iy wings, it is apparently not true for arbitrary non- 
D> 
It must also be noted that these 


latit 
stationary motions. 
s are limited to the lifting problem in which the 


rturbation potential is an odd function of the z co- 


result 
pe 
ordinate. 


(IV) PERTURBATION POTENTIALS FOR VERTICAL 
OSCILLATIONS OF A RECTANGULAR WING 


A flat rectangular wing of zero thickness is shown in 
Fig. 2. Its aspect ratio, A.R. = 6/x,, is assumed to be 
sufficiently large that BA.R. 2 1 (the condition of inde- 
pendent subsonic edges). In the present section it is 
assumed that the wing is undergoing a slight periodic 
vertical motion such that on the wing 

lim (Ob/Oz) = wy exp (if) (7) 


‘ +0 
where wy is a real constant, which may also be expressed 
as 


Mm = - r No (S) 


where Ap is the maximum instantaneous angle of attack 


UA ” wl 
P(x, vy, + 0, ft) = exp (it) 3 > —— 2 = 
; T 0 B*a’ 


By introducing a set of new integration variables de- 
fined by 


o = (vl’/ Ba") (x — &)\ 


: a (10) 
a cos 6 = (vyl’/Ba*) (y — n) J 


the integration with respect to @ can be carried out, and 
Eq. (9) becomes 


P(x, y, + 0, 2) = 


ApBa* s . ‘ o 
exp (zt) exp(—ta)Jo do (il) 
v 0 M 


k = (M*y/B?U)x (12) 


where 


and where J)(o .\/) is a Bessel function of order zero 
and argument a \/. The quantity & will be referred 
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Fic. 2. Rectangular wing in a supersonic flow 


of the wing. The purely supersonic region ODO’ 
and the mixed supersonic regions ODEA, O’DFB, and 
DEF will be designated for simplicity as regions IT, I, 
III, and IV, respectively. 

By application of Eqs. (3), (6), (7), and (8), it is im- 
mediately found, for a point P in region II, that 


(1/A(x — © eos }(v/ Ba) [(x — £)? — By — n)7] 74 , . 
e\9 9 oat adn (\) 
1/B)(x — &) [(x — §)? — By — n)4]” 


s 


to as the reduced frequency.¢ It is obvious in Eq. 
(11) that the potential in region II is independent of the 
spanwise coordinate. This same expression has been 
previously obtained in several investigations of two- 
dimensional oscillating wings.*® 

An alternative expression for this same velocity po- 
tential in region II may be obtained by integration of 
Eq. (9) in characteristic coordinates defined by 


u (M/28) (&§ — Bn) 
v = (M/28) (£ + Bn) J 


(13) 


The characteristic coordinates of point P will be de- 
noted by (#,, v,). In characteristic coordinates, Eq. 
(9) becomes 

+ Other definitions of the reduced frequency, for example, 
yx /U and yMx/p?U, have also been used. The present definition 


has certain advantages in numerical calculations 


i A UAgo - ma du 
(Up, +0, t) = exp (1vt) - 
mM vp (Mp — U) “* exp [(tv/ Ba) (up, — u)| 


By introducing new integration variables 


f° cos |(2v/MBa) [(up — u) (vy — v)] "4 


Mu (v» —v) “exp [tv/Ba) (v, — 2v)| 
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w= [(u — u)/(up + 
P 


it is possible, after some algebra, to reduce Eq. (14) to 


1 
P(t, Up, + 0, t) = (AoBa?/v) exp (ivt) (2k'/*/) f t(u, M, k)du 
0 


— 2k 2 
f(u, M, k) = exp ( ls . ) f 
M? ki 


The definite integral in Eq. (17) is a Fresnel integral and is tabulated. 
Since Eqs. (11) and (16) must represent the same function, two integral representations for the function 


. 
/ exp(— ta) Jo(a/M)do 
0 


where 


T(k, M) = 
= (2k 


have been obtained. t 
parts of T to eight decimal places for0 < k < 5and1 < 


turbance potential in characteristic coordinates. 


AERONAUTICAL 


: vp)]'” l 
(v/Ba) “*[(v, — v) 


This function has been tabulated by Schwarz.° 
M ¢ 10. 

For the point Q on the top surface in the mixed supersonic region I, it is most convenient to evaluate the dis. 
If the characteristic coordinates of Q are (ue, ve, + 0), the gen- 
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*— (1/M) (up — u)'| $ (15) 





16) 


ki /2 [(p/M) + (1 — p?)!/2] 
exp (—iP*)dP 17) 


2 [(p/M) — (1 — p?)! 7/2] 


1 
m / flu, M, k)du (18) 
0 


Schwarz computed the real and imaginary 








eral theory states that ' 
( + 0,1) = —* exp (int) i - 
Ua, Ug, ’ = exp (tv; WA ae 
aM v9 (Ug — u)”* exp [(tv/Ba) (ug — u)] 
| cos | (2v/MBa) [(ue — u) (vg —v)|*} i» (1 
FS dv (19) 
-u (vg — v) exp [(iv/Ba) (ve — v)] 
This integral may be reduced by the same process used with Eq. (14) to ' 
(Bly|/x)'/? 
P(g, Ya, + O, t) = (AoBa?/v) exp (art) (2k “*/m) J f(u, M, Rk) du (20 
0 
+ This identity may also be demonstrated by other methods. 
i 
—~— 
I 
It may be noted that Eq. (20) differs from Eq. (16) only (V) Lirt AND PITCHING MOMENT FOR VERTICAL 
in the limit of integration. The |y| is used, since for the OSCILLATIONS OF A RECTANGULAR WING 
| 


particular coordinate system used y < 0 on the wing. 
It may also be noted that on the Mach line separating 
regions I and II Bly = x; therefore, the velocity po- 
tential is continuous (as it must be) across this 
line. It is of interest to that the 
conical flow variable, 8)y /x, appears explicitly in Eq. 
(20). 

From Fig. 2, 
ditions that identical flow conditions exist in regions | 
and III. It is also possible to determine the velocity 
potential @ at a point in region IV by the superposition 


formula 


note typical 


it is easily seen from symmetry con- 


© = P, + & — (21) 
where ®, is the velocity potential [Eq. (20)] for a right- 
hand wing tip, @ is the corresponding potential for a 
left-hand wing tip, and ; is the two-dimensional po- 


tential [Eq. (11)]. 


In the linearized theory, the disturbance pressur 


for an oscillating wing is given by 
ah iv) (2 


(ue + 2”) a (v= 
Ox or) P Ov 


For the lift problem, which is being considered, tht 
perturbation potential, and thus the pressure disturb 
ance, is an odd function of the z coordinate. Thus 
the dimensionless lift coefficient for any area S 1s 


Ap = 


92 


C, = (4/U?S) f S [U(O’/dx) + ive] dS (23 


The corresponding pitching moment (a stalling momett 
is considered positive) is 


Cu = —(4 U?Sx,) re g [U(ob ox) + iv] x dS (94 


where x, (the wing chord) is the extra length used t 
make the moment dimensionless. 
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MOTIONS OF A WING 


In the present example, it is convenient to calculate 


me 
wing as shown in Fig. 3. 
dimensional center section; Region II’ is the portion 
of the tip in which the two dimensional perturbation 
potential exists; and Region I is the tip region in which 
If it is desired to compute pressure 


Region II is the purely two- 


Eq. (20) applies. 
distribution, it is necessary to evaluate the integrals 
defined by Eqs. (11) and (20). For the overall co- 
efficients, it is not necessary to evaluate Eq. (20) ex- 
plicitly. 

For a strip of unit width in region II of Fig. 3, the 
mean lift coefficient (the two-dimensional result) is ob- 
tained by substituting Eq. (11) in Eq. (23). It is thus 
seen that 


ApBa? ae Bot a er 
~ exp (ivt) — EG + | 1 ax | 
y Us, l 0 


(25) 


where &, is the value of k [Eq. (12)] at the trailing edge 
The quasi-steady lift coefficient for this case 


x = Xe 
iS 
C,,* = (4A0/8) exp (art) (26) 
Therefore, Eq. (25) can be written 
C, l 13" 1 
-, = T(k.) + A(t.) (27) 
Ci," if Mw? ’ 
where 
(28) 


“hy 
A(k,.) = / 1(k)dk 


The corresponding quasi-steady pitching moment co- 
efficient is 
Cy,* = —(2Ao/8) exp (irt) (29) 


The substitution of Eqs. (11) and (29) in Eq. (24) yields 


Ci, 2 137 
ar as k.T(k,) — A(k- 3(k, 3 
Cu,* al (ke) — A(ke) + 47, Beh | (30) 


where 


ke 

B(k,) = / kT(k) dk (31) 

Jo 
For region I of Fig. 3, the perturbation potential of 
Eq. (20) must be used to evaluate the mean coeffi- 
cients. Since the spanwise coordinate occurs only in 
the integration limit of Eq. (20), this spanwise integra- 
tion may be carried out easily. The final result is 
— SAg ar ’ iB? ; 
L Bk? xp (at) | RtT(k.) — A(R.) + ye bee) 
(32) 


where 


7 


*1 
H(k) = (2k'*/x) / (1 — w?)f(u, VW, k)dp (33) 


AT SUPERSONIC 


oe 
~ 


SPEED 

















7 \ 4 


Fic. 3. Rectangular wing subdivisions for the computation of 
mean aerodynamic coefficients 


“ 
C(R,) -f kH(k)dk (34) 


The evaluation of the function //(&) will be carried out 
later. The quasi-steady lift coefficient for region I 
was shown by Busemann’s conical flow theory’ to be 
half that of region II; therefore, 


C,,* = (24Ao/8) exp (ivt) (35 
Finally, 


Cy, = 4 


3? 
= | bate — A(k.) + : = ck | (36) 
C.,* k. a, oe 


The corresponding analysis of the pitching moment 


coefficient shows that 


Cy, _ Of... is | 
= — |b H(k.) — Bk.) — Clk D(k. 
Cu” at ) Tape 


(37) 
where 
“he 
D(k,) = / k?H(k)dk (38) 
0 
and the quasi-steady moment coefficient is 
Cu,* = —(2/3)C,,* = —(4A0/38) exp (art) (39) 


For region II’, the spanwise integration may again be 
easily carried out, and it is thus found that the mean 


coefficients for this region are 


Cr, = 2(Cz, + Cu) (40) 
Cu, = 2Cu, + Cw (41) 
where 
Fal . 4+ eon 13° : 
Cu = —Cy,* 53 | ker Uk, — 2B(k.) + ue Fike 
(42) 
(4:33) 


*ke 
E(k.) - | k?T(k)dk 


For the wing as a whole, the lift and moment coeffi- 
cients may be obtained by weighting the coefficients for 


regions I, II, and II’ in proportion to their areas. Thus, 
Cin = Cr, + (Cr; + 2 Cu) (1/BA.R.) (44) 
Crp = Cure + (Cur, + Cu) (1/BA.R.) (45) 


These formulas were derived with the tacit assump- 
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tion that the tip regions do not overlap (i.e., that 


BA.R. 2 2); however, they are valid if BA.R. 2 1. 


NUMERICAL EVALUATION OF LIFT AND MOMENT 
FUNCTIONS 


(VI) 


In order to discuss the results of this analysis, it is 
necessary to evaluate the various definite integrals de- 
fined by Eqs. (18), (28), (31), (33), (84), (38), and (43). 
These integrals may all be expressed in terms of 7(k) 
and of other well-tabulated functions (Bessel func- 
tions and trigonometric functions); however, it is 
probably more convenient in many cases to use direct 
numerical evaluations of the integrals rather than the 
reduction formula. 


M? by (— ik 
< *X — 1 
32 \ exp ) 


Similarly, 


Pe ese: M? j : k , k 
B(k) = 3 k?T(k) + v7 yF exp (— tk) uw \ + (1 
3iM! § 
284 | 
and 
E(k) l psrpy 4 ep ( | F J ()+e 
i(k) = - k 2 exp (—1 J 2 
lle 3B2. uw \u 


M* { 
384 
. MS { 
36 ) 
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The reduction formula for these integrals can be com 
puted in two groups: first the formula for A(k), Bik 
and E(k); second, the formula for //(k), C(g). and 
D(k). To reduce Eq. (28), first integrate by parts 
Thus, 


a 
A(k) = kRT(k) — | a exp (—t0) Jo(o/ M) da (4g 
0 


Next, replace the Bessel function by its differenti, 
equation, 


wg) wa leda(i)] 


Two additional integrations by parts yield the resy 


l me, 
) — 7 n( )| — iT(k) ¢ 1 


; k eae 
=) 7 (* micelle 
+ in(*) —1T(k) 


d 


do 7 


oO 


M 


k 
M 


M 


wi) 
—_ it) Jo ei 
M 
= ik - 
M M 


al 
yy" 


: (j,) | , ) 
Ji — "T(k)¢ + 
M f 
’ k ee 
) -+- ids ( )| _ iT (k) 


M 


kh 
M 


Reduction formulas of this type have been used by Borbely}* and by Miles.’ 








i 


J 


| 


I 


The reduction formulas for the integrals involving H(k) form the second group. By direct differentiatiw! 


from the second form of Eq. (18), it can be shown that 


of ip? 
oR?) — 


MP 
If this result is substituted in Eq. (33), it is seen that 


d 
dk 


M? 
H(k) = T(k) — Oba? \F exp (—1k) iso ( 


ft °F (u, M, k)du + — ($f E (= 
ro » i, Xx — 12) St 
Tila ne M 


kb 
M 


) 7 : ‘ ( 
)| - iT(h)| 


) 


b 
M 


k 
M 


) a ‘ 2 ( 


{} There is apparently a typographical error in Borbely’s formula. 


This result may be combined with Eq. (48). Thus, 


2kH(k) = A(k) + RT(k) 


(53) 


A substitution of Eq. (53) in Eq. (34) and an integra- 
tion by parts show that 


C(k) = (1/2)kA(kR) (54) 


Similarly, 


D(k) = (1/4) [R?A(Rk) + E(R)] (55) 


These reduction formulas make it possible to co! 
pute all the required integrals from Schwarz’ tabulat 
values of 7(k).° For the numerical results present 
here, a somewhat simpler procedure was used. 7 
function //(k) was computed from 7(k) by Eq. (9: 
Then the related integrals, A(k), B(k) .. . E(k), wet 
computed by numerical integration. The 
formulas were used to spot check the accuracy of t 


reducti 


integration. For certain cases involving small valt 
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imber of significant figures in Schwarz’ tables (VII) PriTcHING OSCILLATIONS OF A RECTANGULAR 


of k, the - . - y - 
inadequate, and supplementary calculations of the WING 
was Iilte ie 
function T(k) by series expansion methods were re- It is well known that the torsional oscillations of a 
quired. In order to illustrate the nature of the re- wing can be decomposed into a vertical oscillation 
ulations were carried through for two Mach (which was considered in the previous sections) and a 


sults, cale 


M = 2and M = 10/7. Only the final re- second oscillation for which the vertical velocity on the 


Numbers : — 
ults required for discussion of the lift and moment char- _ plane of the wing is 

. 7c ¥ > Te ¢ "1 y € > 2S > > > » . - F oe 
acteristics of the rectangular wing are presented here lim (O/dz) = —AUx exp (it) (56) 
in Tables 1-4. More complete tabular data and vari- s= +0 
ous extensions, both theoretical and numerical, of The constant Ay is not necessarily related to the corre- 


arz’ treatment of the function 7(k) are presented sponding constant of Eqs. (7) and (8); however, it is 


Schw 
convenient to use the same symbol. 


in reference 11. 


For a point P in region II (Fig. 2), the perturbation potential is 


ay 


x 
@ = (UAo/m) exp (it) / tg(x — &)dE (57) 
0 


where 


’ iwU(x — &) *y + C/A — 8 eos | (v/ Ba) [(x— £)? — B2(y — n)*] 74 7 
g(x — ~) = exp| — a = eas dn (58) 
82a? ee [(x — £2 — BxXy — 9) 


¢ 


From Eq. (57), — By; however, the upper limit is 7 = —y — (1/8) X 
: (x — £)ifx — — 2 — By. With this modification, the 


- = Uke exp (it) | x«(0) + | . 2 28 ae| (59) argument leading to Eq. (61) is still valid; therefore, 
Ox ” 0 Ox Eq. (61) can be used to compute the perturbation po- 
tential at points in region I (Fig. 2). 

Since the perturbation pressure is defined by Eq. 
(22), the perturbation pressures for the two types of 


Since Og Ox = —(Og/OE), Eq. (59) may be integrated by 


parts, and it is thus seen that 


Ob Uy i er oscillations are also related by a similar formula 
= exp (vt) g(x — &)dé (60) ; 
Ov T 0 7. 
ae ; it : Apr = / Apdx (63) 
rhe right-hand side of Eq. (60) is exactly the same as 0 
that of Eq. (9); consequently, = - ia , ; 
| ? The lift and moment coefficients for the wing are ob- 
" tained by integrating the perturbation pressures over 
*r = a Peds (61) the wing. By introducing Eq. (63) and suitably in- 


verting the order of integration, it is thus possible to 


where &; is the potential defined by Eq. (57) and ®, show that 


is the potential defined by Eq. (9). By substituting “e 
Eq. (11) in Eq. (61), it is readily seen that Cup = (1/2) / x C,dx 
0 a ’ 
sae ™ . : a (64) 
Py = (Aga*B*/ Mv?) exp (ivt)A(k) (62) , a 
Curr = (1/%,*) SC, Ox 

rhis result may also be verified by integration in char- . 
acteristic coordinates. where, again, the subscripts 7 and w have been used 


The formula of Eq. (61) can also be shown to hold for to denote the cases corresponding to Eq. (56) and to 
apont Qin region I. For this case Eq. (57) still holds, Eq. (7), respectively. In Eq. (64), C;,, is the mean 
except that the upper limit in the integral of Eq. (58) lift coefficient for a vertical oscillation of a wing of 


ismodified. The upper limit is as stated if0 < x — € < chord x where 0 < x S x,. 


The lift coefficient C,,.is obtained by substituting Eq. (44) in Eq. (66). Thus, 


C, l 13? } 
“= ) A(k,) + = [k.A(k.) — B(k,)] : + 


rC,.* b2) 


2 iB? 
cea co — Bik.) + = [k.C(k.) — k-B(k.) + E(k.) — Dk (65) 


~9r 
WOO 
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The corresponding pitching moment coefficient is 


Cur l J iB? 
a = < Bik, k2A(k.) — E(k.)] ¢ 
x.C1,* k.3 ( ) = TE [ (R.) ( iF re 
2 J 
kABA.R. | 
where 


“he 
F(k,) J k®T(k)dk 

0 
The function F(R.) may be evaluated by numerical 
integration or by reduction formula similar to that 
used for B(k,) or E(R,). 

(VIII) Discussion 

The fundamental results of the present analysis of 
the vertical oscillations of a rectangular wing, Tables 
1 and 2, are presented graphically in vector diagrams 
in Figs. 4 and 5. The most interesting feature is the 
change in sign of the imaginary part of C,,/C,,* and 
Cy,/Cu,* for M = 2 from that for M = 10/7 and from 
the two-dimensional results. The simplest way of 
demonstrating this effect is to obtain the Taylor series 
expansion (in powers of k,) of the various functions. 
The first two terms of the expansions of the functions 
T(k), A(k),..., E(R) can easily be obtained from their 


integral representations. It is thus found that 


C,,/Cri,* = 1 — (tk,/2M?) + O(k,?) 
Cu,/Cu,* = 1 — (2tk,/3M?) + O(R,?) 
ce. _(M — _ . 

F = 1 ke + O(k-? 7 
c.° te ( 32 al (67) 
Cy _[3(M? — 3) 

= | k. + O(R,?) 

t.* * i Sir | sites 


In each case it is apparent that the quasi-steady limit 


is obtained in the limit as k, ~ 0. It is also seen that 
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Values of the Lift Coefficients Cry and C,, 


k (Cx. / Ci”) 
) 1.00 

0.4 0.980753 
0.8 0.927111 
..2 0.850253 
1.6 0.765481 
2.0 0. 688280 
0 1.00 

0.4 0.990158 
0.8 0.962459 
1.2 0.922000 
1.6 0.876016 


2.0 0. 832383 


Note: The subscripts r and 7 denote, respectively, the r 


ES 


r (CL, 


0 
—(0 
—(Q. 
—Q. 
—Q. 
—0 


0 
—O 
—(). 
—O0 
—O0 
—(), 


TABLE | 


CEs ( 
M = 10/7 


095103 
173680 
223212 
238042 
220289 
M = 2.00 


048598 
089129 
115157 
123145 
113036 


SEPTEMBER, 


: 1B 
D(ke) — Elk.) + 7 Re A(ke) — 3k2B(R.) + 2F(k)) . 
dy. a 


l r 


Cz. /Gs,” 


—_ 


00 

0.977930 
0.917796 
0.836024 
0.754666 
0.695152 


1.00 

0.995200 
0. 982946 
0.969264 
0.956052 
0.972402 


part and the imaginary part of the given function. 


TABLE 2 
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C1,/Ci," 


0 
— (0587 
— 0.09409 
— 0.0927] 
—U 04620 
+) O88: 


0 

0.0347] 
0. 07698 
0), 13235 
0, 20237 


0). 2834 
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Values of the Moment Coefficients Cy, and Cy, 
k (Cu,/Cu,*)r (Cuy/Cu,*)s (Cu,/Cu,*)s (Cu,/Cu,* 
M = 10/7 
0 1.00 0 1.00 () } 
0.4 0.971190 — 0.126030 0.973515 —().0655 
0.8 0.891572 —(). 225748 0.902019 —(). 10318 
‘2 0.779638 — (0.279664 0.806421 —(). 09317 
1.6 0.660424 —(), 279950 0.714660 —() 02058 
3.0 0. 558689 — (0.231794 0.653070 +-() (788 
- - M = 2.00 ——— - 
0 1.00 0 1.00 () 
0.4 0.985260 —().064396 0.994125 0.03% 
0.8 0.944098 —(. 115986 0.979755 (0) O878% 
12 0.884984 —(0. 144622 0 (0). 1527 
1.6 0.819668 —0.145192 0.9: 0). 22741. 
2.0 0.761014 —0.118418 0.97557 () 33223 
04— ; | 
sa —_—_ M= 10/7 
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TABLE 3 
Lift Coefficient of the Rectangular Wing, Cz, 


SONIC 


SPEED 





- - M=10/7 : 
A.R. = 1/B A.R. = 1 A.R. = 3 AR. =5 
b Cr./Cu.)y (Cx,,/C1,)s (Ci, /C1,)r (Ci,,/C1,)i (Cr,,/C1,)r (Cx, /Ct1,): (Cx,,/Cug)r Cr /Cs 
0) =X) 0 0. 509902 0 0.836634 0 0.901980 0 

- 0) 503422 0.050418 0.513256 0.049420 0.837752 0.016473 0.902651 009884 
(0 : 0) 514336 0. 101284 0.528954 0.099278 0.841318 0.033093 0.904971 0). 019856 
re 0) 534715 0.152248 0.543930 0.149233 0.847977 0.049744 0. 908786 0.029847 

f 0) 567663 0.201032 0.576225 0.197051 0.858742 0.065684 0.915245 0.039410 
. > 0.615924 0). 241677 0.623530 0.236891 0.874510 0.078964 0.924706 0.047378 
: M = 2.00— 
() 0). 500000 0 0.711325 0 0.903775 0 0.942265 0 7 
04 () 504634 0.058250 0.714000 0.033631 0.904667 0.011210 0.942800 0. 006726 
08 (). 518981 0.115958 0.722284 0.066948 0.907428 0.022316 0.944457 0.013390 

9 (0). 544384 0.171726 0.736950 0.099146 0.912317 0.033049 0.947390 0.019829 

f 0). 578792 0.222041 0.756816 0.128195 0.918939 0.042732 0.951363 0.025639 

) 7 ’ ay eee ome K“A\-oO —— me Owe 
» () 0.634156 0). 262832 0.788780 0.151746 0). 929593 0.050582 0.957756 0.030349 

TABLE 4 
Moment Coefficient of the Rectangular Wing, Cy, 
AR. = 1/8 A.R. = 1 A.R. = 3 AR. = 5 
b Cu /Cm (Cu /Cu (Cu, Cu, )s ‘Cu, Cm, ): (Cm, C,)r (Cm, Cm,)i (Cy /Cu Cu /Cu 
- M = 10 ‘ = a 

( 0.333333 0 0.346536 0 0.782179 0 0. 869307 0 

} 0.334598 0.062541 0.347776 0.061302 0.782592 0.020434 0.869555 0.012260 
08 0.339791 0). 128996 0.352863 0.126441 0.784288 0.042147 0.870573 025288 

9 0.341073 0.199478 0.354122 0.195528 0.784707 0.065176 0. 870824 0.039106 
1.6 0). 381829 0.289449 0.394071 0). 283717 0.798024 0.094572 0.878814 0.056743 
20) 0.445849 0.382746 0.456823 0.375166 0.818941 0.125055 0.891365 0.075033 

M = 2.00 - 

() 0). 333333 0 0.615100 0 0.871700 0 0.923020 0 
04 0.338132 0.081429 0.617871 0.047013 0.872624 0.015671 : 0.009403 
0.8 0.353718 0. 163898 0.626869 0.094627 0.875623 0.031542 0.018925 

2 0.384502 0.249667 0.644642 0.144145 0.881547 0.048048 0.928928 0. 028829 
1.6 0.481251 0.327477 0.671633 0. 189069 0.890544 0.063023 0.934327 0.037814 
20) 0.523049 0.410011 0.724632 0.236720 0.908211 0.078907 0.944926 0.047344 


the two-dimensional solutions show a lagging phase (the 
imaginary part is negative) for small values of the re- 
duced frequency for all Mach Numbers. On the other 
hand, the solutions for the tip region show a phase 
lag for 1 < JJ < 3 * and a phase lead for J > 3 

For the wing as a whole, it is probably best to consider 
the ratio of the three-dimensional and the two-dimen- 


sional results. From Eqs. (44) and (45), 


Cue _ 1 S _ | | 

Cz Cio _ BA.R. (68) 
“Mu = | , (‘ Vv + “*) l 

Cy Cw. BA.R. 


The computed values of these functions are listed in 
ables 3 and 4 and are plotted in vector diagrams in 
Figs. 6-9. These results show that in every case the 
effect of the tip is to shift the phase of the lift and 
moment in the positive (leading) direction. By com 
parison with the two-dimensional results of Figs. 4 and 
), it is apparent that in many cases the effect of the tip 
Is strong enough so that the sign of the phase is re 
versed. Since the phase angles are extremely im 
portant for flutter calculations, it is apparent that 
important finite span effects are involved. 


A second method of examining the overall effects 
would be to consider the ratio of the three-dimensional 
results to the two-dimensional quasi-steady results. 
For this comparison, the Taylor series expansions are : 
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stability calculations, the results given in Eq. 
should be sufficiently accurate representations of th | 
nonstationary effects. , 
Some of the results of the pitching oscillation ana}, 
can be seen from the Taylor series expansions of p 





(65) and (66). 


ti 7 (1 | ) 
x.C,.* 9 3BA.R. 


ik, ( M? + ) 
(1 + O(R2 
6M? 


These expansions are 


4BA.R. 
— Cur (i _  B ‘- 
wie.” 3 SBA.R., 
tk, 1M? 4 
= (: —— : ) + O(R, 
SA? 15 BA.R. 


It may be noted that 


x-Cy,* = (4x-A0/B) exp (it) 


which is the two-dimensional lift coefficient based 
[see Eq. (56)] the instantaneous angle of attack att 


trailing edge. It may be noted that these solutic: 





like those for the vertical oscillations, show phase 1 
versals for finite aspect ratio if M/ is large enough 
One of the most interesting results of the two-din 


sional theory was the one degree of freedom torsio: 


It is of 
terest to note the modifications to this theory int 
If the wing 
executing a torsional oscillation of amplitude A, ab 


flutter noted by Garrick and Rubinow.* 
duced by the finite aspect ratio effects. 


the chordwise point x9, the position of the points ont 
wing surface are given by 


n(x, t) = —Ao(x — Xo) exp (at 


The required vertical velocity on the wing is then 


. Of 4 1vXy w\y | 
lim = =f exp (zvt) | Ao (: — : t \ 
5 0 OZ [ l J 


~ 
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Ce (: l ) 
Sky 2BA.R. 
tk, M? + 1 , 
rae + ((k,7) 
2M? 3BA.R. ' 
(69) 


City _ (1 2 ) 
Cu.* 3BA.R. 
21k. 3M? 3 
(: — Z ) + 0(k,2) 


31 SBA.R. 
From this result it can be seen that for any aspect ratio 
there will be a range of Mach Numbers for which the 
three-dimensional wing will have, at least for small 
reduced frequencies, a phase angle that is opposite to 
that of the two-dimensional theory—i.e., which leads 
the quasi-steady value. For many purposes, such as 
the calculation of aerodynamic damping for dynamic 
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TABLE 5 
Lift and Pitching Moment Coefficients for a Pitching Oscillation 
k Ly Li, Le, Le, M,, M, Mo, Wo. 
Z M = 10/7 
() 0). 500000 0 —(). 166667 0 —0. 333333 0 0.125000 0 
04 (). 495158 — (0.032088 —(. 163050 0.024875 — 0.329453 0.024016 0.121954 —(). 019532 
V8 (). 481325 — (0.060806 —0. 152922 0.046228 — 0.318449 0.045234 0.113560 —(). 036802 
. (0) 460434 — (0.083380 —(0.137658 0.063408 —(0.3019383 0.061367 0.100958 —(Q) 050122 
& 0.435269 — (0.098067 —(0.119478 0.074694 — 0.282282 0.071029 0. 086090 —0. 058468 
90) (). 408935 —(). 104892 —0. 100673 0.079807 — (0.262107 0.073932 0.070890 —() 061683 
i M = 2.0 
() 0. 500000 0 —(). 166667 0 — 0.333333 0 0. 125000 0 
4 0.497528 —(.016400 —0.164225 0.020326 —(0.331375 —(.011609 0.123750 0. 016641 
V8 () 490410 —(0.031136 —(.157052 0.039631 —(). 325689 — (0.023174 0.116987 0.031587 
» =). 4795383 —(). 042846 —(). 146018 0.055823 — 0.317073 — 0.031560 0.107856 0.044280 
16 (). 466182 —(). 050549 — (0.132322 0.070879 — 0.306617 — (0.036649 0.098158 0.056414 
» () 0. 451876 —(). 053827 —(0.117606 0.075879 —(). 295592 —(0) 038142 0.084163 0.058911 
The complete solution is thus obtained by adding to- Tene 6 
vether the vertical oscillation of Eq. (69) } with Ao re- Lift and Pitching Moment Coefficients for a Vertical Oscillation 
placed by Ao[1 — (ivxo/U)}} and the torsional oscilla- k L L M V 
tion of Eqs. (70) and (71) [with Ay replaced by (iv, U’) X M =10/7 
Ayl. The pitching moment coefficient with respect to 0 —0.500 0 —(0).666667 0 
tae AG C sa teas 0.4 —(). 482225 0.096675 —(). 638350 0. 144600 
the axis of rotation © M;, 1s thus 0.8 —(). 432674 0.178253 —(0). 559504 0. 264050 
; 1.2 —(). 361626 0. 233307 —(). 447936 0.339798 
Cy, l l Xo I 1.6 — (0. 283091 0. 256802 — 0.327226 0.364216 
ce —\o~ 30AR + - _- 2BA.R 2.0 —(. 211113 0.250949  —0.220880 0.342284 
L yA OPrr. . “We wll. ° Vv - 90) 
1k, | (? — I ) 0) —0. 500 0 —0. 666667 0 
2 “ cat is 0.4 —(). 487660 0.081751 —(). 646870 0. 122850 
M Ki IBA.R. 0.8 —().452625 0.154479 —().591146 (0). 229696 
Xo (2M? — 3 9M? — 4 ie —(). 400302 0.210800 —(). 508568 0. 309954 
— — 1.6 —(0. 341642 0.246380 —0.418640 0.351002 
+. 2 6BA.R. 2.0 —(). 274813 0. 260130 —(0.314414 0. 368502 


(“) (: l ) e| + Ok) (74 
ma 2BA.R. os os 


The terms that are independent of k, give the quasi- 
steady static stability moment. The next term gives 
the damping in pitch, provided the reduced frequency 
is small enough so that the higher order terms in k, 
may be neglected. 

The motion is damped if the quantity multiplying 
tk, is negative; the negatively damped, one degree of 
freedom flutter condition occurs if the quantity multi- 
plying 7k, is positive. For infinite aspect ratio, Eq. 
(74) shows results identical to those of reference 6. 
For example, it is seen that an instability is possible, 
with suitable values of xo, if 0 << B2?< 1.5. If B? = 1.5, 
the motion is neutrally damped if xo/x, = 1/3. If 6? = 
|, the motion is negatively damped if 0 < xo/x, < 1/2. 
For finite aspect ratios the region of instability is con- 
siderably reduced. If 8 = 1 (M = 1.414), instability 
is possible if A.R. > 3.414; if 8 = 0.707 (M = 1.224), 
instability is possible if A.R. > 1.99. It thus appears 
that, for the small aspect ratios that are normally 
used, the range of Mach Numbers for which this type of 
instability is theoretically possible is in the upper 
transonic range where the linear theory is not reliable- 
however, the theory does indicate that the usual method 
of estimating aerodynamic damping may be seriously 
in error 

For finite values of the reduced frequency, the calcu- 
lated values of the lift and moment coefficients for the 


pitching oscillation can best be tabulated by writing 


Eqs. (65) and (66) in the form 


= Li(k,) + [L2(k-) BA.R.] 
= M,(R,) + [ Mo(R,) BA.R.| § 


i ss 
C1,/x clr, 
(75) 


. we 
—( MT x Lo 


The real and imaginary parts of the functions L,(&,), 
Lo(k.), My(k-), and M.(k,) are given in Table 5. In 
order to simplify the computation of vertical oscilla 
tions, Eqs. (44) and (45) can be rewritten in the corre- 


sponding form 


(CL, C.,*) + [L3(R,) BA.R.| ' 
*) 4 [My(k,)/BA.R.| § 


Cs. tg _ nam 
. + * . . (460) 

tate Cx = (Cus, Cu 
The two-dimensional values (8A.R.— ©) are tabulated 
in Tables 1 and 2. The functions L3(k,) and /3(k,) are 


tabulated in Table 6. 
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A Recurrence Matrix Solution for the 
Dynamic Response of Elastic Aircraft 


JOHN C. HOUBOLT* 
Langley Aeronautical Laboratory, N.A.C.A. 


SUMMARY 


A systematic procedure is developed for the calculation of the 
structural response of an airplane subject to dynamic loads. 
Particular attention is given the problem of determining the 
stresses developed flight through gusts. Difference 
equivalents for derivatives and matrix notation are used to de- 


due to 


velop a recurrence relation that permits step-by-step calculation 
of the response and of the loads that occur on the structure. 
The chief feature of this recurrence approach is that the general- 
ity and physical aspects of the basic equilibrium relations of the 
problem are preserved without loss of ease in application. 

The use of difference equivalents for derivatives in the solu- 
tion of dynamic problems is first illustrated by means of a simple 
damped oscillator example, and the application to the flexible 
aircraft structure is then made. For brevity, the case of wing 
bending and vertical motion of the airplane is treated, although 
the method may be readily extended to take into account also 
wing torsional deformations, pitching motion of the airplane, 
fuselage deflections, and tail forces of known character. Either 
a sharp-edge gust or a gust of arbitrary shape in the spanwise or 
flight directions may be treated. 

Some results obtained by application of the recurrence ma- 
trix relation are presented, and the advantages of this method 
over other methods of evaluating the dynamic response of an 
aircraft are discussed. 


SYMBOLS 


E = Young’s modulus of elasticity 

I = bending moment of inertia 

m = mass of beam included in the interval / or concen 
trated mass in spring oscillator 

m = mass m including apparent mass effect, [m + 
(wplc?/4) | 

a = distance between leading edge of wing and elastic 
axis 

r = chord of wing 

Co = chord at wing midspan 

b = semispan of wing 

l = length of section associated with a spanwise station 

Ni = dimensionless interval between the 7 — 1 and ith 
station (A;) is actual length) 

A = aspect ratio of wing 

Ww = deflection of elastic axis of wing, positive upward, 
or deflection of mass oscillator 

g = ratio of dynamic deflection to static deflection 

¢ = angle of twist, positive in stalling direction 

tir = time, zero at beginning of gust penetration 

€ = time interval 

§ = distance traveled by wing in half-chords [(2U’/co)f, 
where midspan chord ¢ is used as reference 
chord | 

r = exponential coefficient in the @ function associated 


with variable s 
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¥ = exponential coefficient in the # function associate 
with time ¢, [y = (2U/c)A] 

n = integer 0, 1, 2, 3, etc., used as subscripts to desi 
nate number of time intervals passed 

1(t) = unit-step function 


F = suddenly applied force or portion of the forces dy 
to the wake 
mass density of air 








p = 

U = forward velocity of flight 

M = Mach Number 

i = vertical velocity of gust 

m4 = assumed overall compressibility and aspect-ratj 
correction, [4/(2 +A V1 — M?)] 

B = forward speed factor for wing (m4 7p U) or coefficien: 
of damping for spring oscillator 

1 — # = function that gives growth of lift on airfoil following 
sudden change in angle of attack 

y¥ = function that gives growth of lift on rigid airfo 
entering sharp-edge gust 

p = normal load acting at a station 

a = aerodynamic lift over interval / on wing due to gust 


square matrix 
= column matrix 
Dots are used to indicate derivatives with respect to time; for 
example, Ow/dt = ¥ or Ow/Or = w. 
Station locations are denoted by use of parenthetical numbers 
(except for \;, where 7 denotes the station). 


INTRODUCTION 


y ] NE POSSIBILITIES THAT critical loading conditions 
may be encountered by an airplane in flight through 
necessitated the development of 


gusts has better 


methods for predicting the stresses that are induced 





In the case of flexible aircraft, accurate prediction oi | 


stresses cannot be obtained if the interaction between 
aerodynamic loads and structural deformations is not 


considered. 


The present paper presents a method for determining ! 


the dynamic response of aircraft in gusts in which 


this interaction is considered. An approach is em 
ployed which is a departure from the usual modal type 
of solution. The time the integro 


differential equations of motion of the airplane are re 


derivatives in 


; : 
placed by appropriate difference expressions, and ust 


is made of matrix notation to express conveniently the 
conditions of equilibrium at a number of points along 
the wing span. The result is a rather simple system 
atic procedure, which is, nevertheless, complete and 
general in form. Freedoms of the airplane in vertical 
motion and pitch, flexibility of the wing in bending 
and torsion, flexibility of the fuselage, and tail forces 


may all be included without difficulty. For brevity ™ 
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DYNAMIC RESPONSE 
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(0) Damped oscillator and suddenly applied force. 
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(b) Response. 


Figure |.-Exact and difference equation solutions for 
response of damped oscillator. 





presentation, however, the present analysis is restricted 
and wing bending. 

a gust with any gradient in the 
along the span may be treated. 


to vertical translation 

With this method, 
direction of flight or 
The method is based on aerodynamic strip theory, 
but approximate overall compressibility and aspect 
ratio corrections may be included without difficulty if 
correction is indicated. 


desired. One such overall 


In the first part of the paper the method of using 
difference equations in the solution of dynamic prob- 
lems is illustrated by an example in which the transient 
response of a simple oscillator is determined. The 
analysis for the determination of the response of an 
airplane in a gust is then given. 

Supplementary definitions and derivations are pre- 
sented in appendixes. Appendix A summarizes the 
various matrix coefficients and matrices that are used 
in the analysis, Appendix B gives the derivation of the 
flexibility matrix, and Appendix C gives a derivation 
of a recurrence equation for evaluating the form of 
Duhamel's integral that is encountered in the analysis. 


TRANSIENT RESPONSE OF A SIMPLE DAMPED OSCILLATOR 


In order to illustrate the use of difference equations 
and to test the accuracy of the procedures that are to 
be used in the solution of the more complicated gust 
problem, the solution of a simple problem having a 
known analytical solution is presented first. The 
problem is to compute the response of the damped os- 
cillator shown in Fig. la to a suddenly applied force. 
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The differential equation of motion of this system due 
to the suddenly applied force is 


mw + Bw + kw = FI(t) (1) 


By use of difference equations, this differential equa- 
tion may be transformed into an equation that involves 
deflection ordinates at several successive values of 
time. 

Probably the most commonly used difference equa- 


tions are the following: 


Wy = (Writ — Wp_-1)/2E (2) 


Wy = (Waa — 2w, + Wr-1)/e (3) 
which give the derivatives at the intermediate of three 
successive ordinates (see Fig. 2a). Although these 
equations are adequate for the oscillator problem of the 
present paper, they are unsuitable for use in the tran- 
sient response determination for beams for reasons that 
will be brought out in a subsequent part of the analysis. 
For application to beams, it is expedient to introduce 
the derivatives at an end ordinate, and to do this a 
fourth ordinate is added (see Fig. 2b) so that the re- 
sultant derivatives have an accuracy roughly com- 
parable to Eqs. (2) and (3). Thus, from a considera- 
tion of a cubic curve that passes through the four suc- 
cessive ordinates, the following difference equations 
for end ordinate derivatives may be obtained: 


W, = (1/6e) (11 w, — 18 Wy_-1 + DwWy-2 — 2Wy-s) 


(5) 


w, = (1/e*) (2w, — SwWy-1 + 4Wyp_2 — Wy_3) 


Of use also are the derivatives at the third of the four 
ordinates given by the following equations: 


w, = (1/G6e) (26,41 + 3w, — GWyr_-1 + Wp-2) (6) 


W, = (1/e7) (Wasi — 2Wy + Wp-r) (7) 


For illustrative purposes, Eqs. (4) and (5) rather than 
Eqs. (2) and (3) will be used in the solution of the 
oscillator problem. 

If the derivatives in Eq. (1) are replaced by the dif- 
ference Eqs. (4) and (5), the following equation is ob- 


tained: : 


llBe k . _  BBe 
i ae +—é€)w, =\9+ Wy—1 
6 m m m 
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(b) Cubic 
Figure 2.- Notation used in derivation of difference equotions. 
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Figure 3.- Division of wing into sections. 





This equation may be said to be the difference equation 
of motion. 
applied force were being considered, the factor F in this 


If the more general case of a variable 


equation would be replaced by F,, the value of the force 
p-esent at the time ¢ = ne. 

If a specific case is now considered in which k/m = 
400, B/2m = 2, ¢« = 0.01, F = 1, and the notation z = 
w/(F/k) (ratio of dynamic deflection to static deflec- 
tion) are used, Eq. (8) becomes 


1.921142,. + 
0.479492,,3 


Zn, = 0.018927 + 2.42272z2,_, 
(9) 


This equation may be regarded as a recurrence for- 
inula; the value z, may be interpreted as the deflection 
to come and may be found easily from the three pre- 
ceding deflections 2,1, 2,2, and 2,3. Then, with 
the newly found value 2, and with 2,_; and 2,_», the 
next deflection can be found and so on. This process 
thus gives a step-by-step derivation of the time history 
of deflection and may be carried out as far as desired. 
Of course, the process must start with known initial 
values of z. These values can be found with the aid of 
the initial conditions of the problem by means of Eqs. 
(6) and (7). 


t = 0(n = 0) they become 


If these equations are taken to apply at 


(1/G6e) (2m, + 3wy — Bw_y + w_s) (10) 


Wo = (1/e7) (wi — wo + w_s) (11) 

For the problem under consideration, the primary 
initial conditions are that at ¢ = 0 the displacement and 
velocity are zero. By use of Eq. (1) or by reasoning 
from Newton's second law, a secondary initial condi- 
tion can be established—that is, the acceleration im- 
mediately following the application of the unit force 


must be 1/m. In equation form, these conditions are 


W = 0, WwW = 1/m 


By substitution of these values into Eqs. (10) and (11) 
and by use of the notation for z, the following relations 
can be found to exist between the ordinates: 


29 = 0.24 — S21, 21 = 0.04 — 2) 


Za = O, 


Substitution of these values into Eq. (9), with m set 


equal to 1, gives an equation from which 2), the de- 


flection at ¢ = ¢, may be evaluated. Three successive 
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deflections can now be established: the deflectioy , 
0, and a fictitious « 
—e given by Eq. (12). In the real 
problem no deflection exists for ¢ less than 0; the a¢. 


: t 
t = e, the zero deflection at t = 


flection for t = 


sumption that a deflection does exist before 1 is zero ic 





simply a means for allowing the recurrence formy 


[Eq. (9)] to apply at the origin as well as at later 


values of time. Furthermore, no violation is made oj 
the conditions under consideration because, mathe 
matically, the response after ¢ = 0 is not influenced }y 
the response that may exist before ¢ = 0 so long as thy 
initial conditions are satisfied. The process just ¢ 
scribed for treating the initial conditions is actually yy 
different from the procedure commonly employed in 
difference equation approaches in which exterior points 
near the region under consideration are written in terms 
of the interior points by means of the boundary cond; 
tions. 

With the initial values of deflection thus established, 
the step-by-step evaluation of succeeding deflections 
proceeds in a straightforward manner—that is, Eq. (9 
is evaluated for m = 2, then for n = The re 
sponse of the oscillator for the physical constants listed 


3, etc. 
previously is given in Fig. lb. The comparison be. 
tween the difference solution and the exact solution js 
seen to be good. 

RECURRENCE EQUATION FOR RESPONSE 0! 
AN AIRPLANE IN A GUST 


MATRIX 


The application of difference equations to the response 
determination of an airplane due to gusts or other 
transient loads may now be shown. The case treated 
is an airplane having an essentially straight wing flying 
through a gust of known structure. For the sake of 
brevity, only wing bending and vertical motion of the 
airplane are considered. The extension is readily made 
to the more general case where pitching motion, wing 
torsion, and fuselage bending are also included. 

Among the assumptions made in the present analysis 
are the following: (1) Engineering beam theory can be 
used to calculate wing bending; (2) aerodynamic strip 
theory applies. effects, however, 
are taken into account by means of approximate over- 


Three-dimensional 


all corrections. 

It is considered convenient to divide the problem 
into two phases, one being the determination of the 
wing deflection under applied loads and the other being 
the determination of the loads from the dynamics of the 
airplane motion. For the determination of the de- 
flections, the wing semispan is considered to be divided 
into a number of sections (not necessarily equal) with 
a station point at the center of each section (see Fig. 
3). In the been 
adopted, although more or less can be used as desired. 
Station 0 is located near where the fuselage intersects 
the wing, and the other stations are located in any con- 


case shown, six stations have 


venient manner along the semispan so as to fall at con- 
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Figure 4.- Forces on an airfoil in transient motion. 








centrated mass locations or at the average of dis- 
tributed masses, station 5 being located nearest the tip. 
The mass within and the air force that develops over 
the section are considered to be concentrated at the 
station point, and the average chord of each section is 
assumed to apply. In this way, the wing semispan is 
considered to be one-half of a free-free beam subjected 
to six load concentrations and, as such, will have a linear 
moment distribution between each station. With the 
further assumption that the |/// variation is linear 
between each station, equations may be developed con- 
veniently by analytical methods which allow the cal- 
culation of the deflection at each station point. These 
equations are developed in Appendix B and may be 
expressed by the following matrix equation 


[A] |w| = |p| (13) 


where [A] is a flexibility matrix involving the // 
values at each station and the station intervals, |w) is 
the matrix that involves the deflection of the elastic 
axis at each of the station points, and |p; is the load 
matrix. This equation thus represents six simultane- 
ous equations allowing the deflection to be obtained 
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Figure 5.- Lift functions for sudden change in angle of attack. 
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directly from the load. The loads are now to be de- 
fined. 

The loads that develop on the wing are defined in 
terms of the inertia forces, the aerodynamic forces due 
to airplane motion, and the aerodynamic forces due to 
Of these, the aerodynamic forces 
are the most complicated. ‘ that an 
airfoil that has freedoms in both translation and rota- 
tion and which is penetrating a gust has acting on it 
aerodynamic forces that may be resolved into the force 
systems shown in Fig. 4. The force L, is the lift force 
developed directly by the gust and corresponds. to the 
gust force that would develop on the airfoil considered 
All other 
These 


gust penetration. 
It can be shown! 


rigid and restrained against vertical motion. 
forces occur as a result of motion of the airfoil. 
forces, as well as the gust forces, are given for an in- 
terval / of the span by the following equations: For the 
forces due to circulation, 


“¢ Ov 
Ly = myrocl f v(t — r)dr (14) 
0 Or 
hs 3 \ 
L, = merocl ff | Ue —wt+ re — -) ‘| x 
0 4 Cc 
[1 — &(t — r)|dr (15) 
Le = (marplc?/4) Ue (16) 


and for the aerodynamic inertia force and moment, 


tplc* r Cc 
4 | -w * (5 


M, = —(nplc! 


L; = (17) 


Lg (18) 
where m, is a factor that can be used to introduce over- 
all compressibility and aspect ratio corrections; in 
this paper the factor is assumed to be given as the prod- 
uct of a compressibility aspect ratio correction defined 
by A’ 
the 


(A’ + 2), where A’ is equal to A V1— M?, 


and Glauert-Prandtl Mach Number correction; 


thus, 


A/(2+ AV1 — M?) (19) 


nN, = 


| — # is the lift function that denotes the growth of 
lift on an airfoil following a sudden change in angle of 
attack; the ® functions may be given by the useful 


approximations (see Fig. 5): 


(1 — 6), = 1 — 0.28%e°°™ (20a) 
(1 — @),-6 = 1 — 0.361e7°*" (20b) 
(1 — ’),-. = 1 — 0.165e~°"” 0.335e7°"""5 ~~ (20c) 


y is the lift function that denotes the growth of lift on 
a rigid airfoil entering a sharp-edged gust and may be 
given by the useful approximations (see Fig. 6): 


Yaes = 1 — 0.679e~°** — 0.22727 °°" (21a) 
Va-s = 1 — 0.448e~°°* — 0.272e7°"% — 
0.193e~*""S (21b) 
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Ws-. = 1 — 0.236e7°°" — 0.513e7°*** — 

O.171e~***  (21c) 
View = 1 — 0.500e~""™”* — 0.500e~° (22) 
Eqs. (20) to (21) are based on equations presented in 


reference 5, whereas Eq. (22) is the y function that is 
suggested for wings of infinite ratio in reference 3. In- 
spection of Eqs. (20) shows that the ® function for as- 
pect ratios 3 and 6 is given by a single exponential term. 
It is probable that the # function for higher aspect 
ratios, say 10 and even 20, may also be given to a suffi- 
cient approximation by a single exponential term, and 
therefore the assumption is made that, in general, ® 
may be represented by an equation of the form 


= ¥ 


—Xs = . 
® aye = aye (23) 


where y (2U/co)X. Interpolation, for example, of 
the curves in Fig. 5 shows that the function for aspect 


ratio 10 might be approximated by the equation 


®,-10 = 0.4le~°* (24) 


It is to be noted, however, that the analysis does not 
necessarily limit ® to a single exponential term. 

Although Eq. (14) gives the gust force for a general 
gust gradient, it is considered worth while to deter- 
mine the response due to a sharp-edge gust, for with 
this response the response due to other types of gust 
may be found directly by superposition. Also, since 
only vertical motion of the airplane is being considered 
in the present analysis, all terms involving the twist 
yg are assumed to be zero. Thus, for the case of the 
sharp-edge gust penetration of an airplane having free- 
doms in vertical motion and wing bending, Eqs. (14) 
to (18) become 


L, = Belop (25) 
ad 

Li = —Bel / w{1 — &(t — t)|dr (26) 

L3; = —(mplc?/4)w (27) 


where § has been introduced as a forward-speed and 
aspect-ratio parameter defined by the equation 


B = m,xpU (28) 


The total load p that acts at a station may now be 
found as the sum of the inertia loads and all the aero- 
dynamic loads that act; thus, 


p= —-moe+1N4+4:+ L, 


The term Z; is associated with the apparent mass effect 
and the inertia term. 
The equation for load becomes then 


p —mi+tl+L, 


may be combined with mass 


= (29) 


where 


m = m+ (plc?/4) 


SCIENCES—SEPTEMBER, 


1950 











on 1 

‘Reference 3 
fe) 4 j 4 1 4 n 4 1 
10) 4 8 12 16 


s, half chords 


Figure6.- Lift functions for wings entering a sharp edge gust, 





Difference Eqs. (4) and (5) may now be used to x 
The term |, 
The $(¢ — ; 
function is associated with the lift forces that are dy 
Without this term the equation woul 


duce the loading to difference form. 
however, requires special consideration. 


to the wake. 
yield the steady-state lift corresponding to the instan 
If Eq. (26) is 
integrated by parts and the initial conditions are used 


taneous value of the vertical velocity. 


such that w and w are 0 at ¢ = 0, the following equatio; 
may be obtained: 


7, 


Li = Belébw — (1 — %)Bclw + Bel / w(t — r)\dr 





where, by Eq. (23), 


Po ay, Dy —(2U Co) AQy 


Cc 


equation may be replaced by the ex- 


In Appendix C, a method is derived whereby the 
integral in this 


pression 
F,, + (1/2) Bcledow,, 
where 
(30) 


F.,, - é as 1 + ZW), 1 


in which g is defined by Eq. (A5) in Appendix A. The 


' 


value of Z; at the mth time interval may therefore be | 


written 


L Bcl [dy + (1 2)eby|w, = (1 = Py) Bclw, Ya F, 


ii = 


With 
be transformed into the form 


the use of difference Eq. (4), this equation may 


Lin = doWy + dW, + down_2 + d3w,_3 + F, (31) 
where the d's are defined by Eq. (A2) in Appendix A. 
With this equation, Eq. (25), and difference Eq. (5), 
the loading Eq. (29) at the mth time interval may be 
written 

P = MW, + Mwy + MWp2 + N3Wr-s + Fy + An 


29) 


where the 7's and f are defined by Eqs. (A3) and (A4) in 
Appendix A and where y, is the value of the y function 
at f ne. This equation thus gives the load in dil- 
ference form and applies at each spanwise station. The 
coefficients y's are seen to involve only the physical 
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s of the airplane structure, the forward speed 


propertie 
aspect ratio parameter given by Eq. (28), certain 


initial values of the unsteady lift function, and the time 
The time interval that is chosen should be 
fairly small in comparison with the natural period of the 


and 


interval e. 
fundamental mode in bending of the wing. An interval 
chosen near one-thirtieth of the estimated natural 
period of vibration of the fundamental mode appears to 
be satisfactory. 

Eq. (32) for loading may now be combined with Eq. 
(13). In order to simplify this process, the abbreviated 
or symbolic form of the matrixes that occur will be used. 
The definitions of these matrixes are given, unless 
otherwise stated, in a complete group in Appendix A. 

Application of Eq. (32) to each of the spanwise sta- 
tions leads to a set of loading equations which may be 
put in the matrix form given by the following equa- 


‘ 


tion: 


Pin = [no] Win i Im] ||» 1 = 5 [n2} lw n—2 + 


[ns] \teln3 + |Fin + \f ln (33) 
where [see Eq. (30)] 
Fin = @ |Flaa t+ ([g] |win—-1 (34) 
Substitution of this equation into Eq. (13) gives 
[A] \w!n = [no] |t\n + [m] ||n—-1 + [ne] w|,2 + 
[ns] |@\n-3 + |Fin + |f itn (35) 
With the use of the notation 
[D] = [A] — [nl] (36) 
and 
Qe = [m] \@\n1 + [me] \win-2 + 
[713] [na + |Fin + \f lve) (37) 
Eq. (35) may be written simply 
[D] |wl,n = IQ, (38) 


Multiplication through by the reciprocal of [D] (the 
Crout method, reference 6, serves as a useful tool to 
make this inversion) gives, finally, the equation 


w\, = [D]-|Q|, 


This equation is the recurrence equation for response 
and gives the displacements that apply at time m in 
terms of the displacements that occurred at several 
preceding values of time [see Eqs. (34) and (37) for the 
definitions of | F| and |Q)]. 

To start the response evaluation, three initial values 
of response have to be known, and these are found from 
The initial response is 


(39) 


the initial conditions as follows. 
assumed to be given in terms of four successive ordi- 
nates, denoted by w_s, w_1, wo, and w,; the wo ordinate is 
considered, as in the case of the damped oscillator, to 
locate the origin of time. The first and second deriva- 
saa at the w) ordinate are given by Eqs. (10) and (11). 
The initial conditions are that the vertical displace- 


ments and vertical velocity are zero. Since the gust 


OF 


BLASTi¢c AtRCRAPrT 545 


force starts from zero, the additional initial condition 
can be established that the acceleration must be zero 
at the start of the response. The ordinate wy and the 
derivatives given by Eqs. (10) and (11) must therefore 
be zero, and thus ordinates w_. and w_; are found to 


to be related to ordinate w, by the following relations: 


(40) 


— SW 


a0) _— 


(41) 


—W 


These relations are general and must apply to the de- 
flection at each spanwise station—that is, the displace- 
ments at ¢ = —2e must be —S8 times the displacement 
at ¢ = e, and the displacement at / = —e must be the 
negative of those at ¢ = e«. Substituting these condi- 
tions in Eq. (35), taking ” as equal to 1, and remember- 
ing the condition that the displacements are zero at t = 
0, the following matrix equation is given: 
[(D] + fm] + Sins] Jeol = [fly = (42) 
The term |F); is zero and therefore does not appear in 
this equation. Solution of the equation gives the val- 
ues of the displacements that occur at ¢ = €(” = 1). 
The three sets of initial displacements required to pro- 
ceed with Eq. (39) are thus known: The set of deflec- 
tions found at ¢ = e, the zero set at ¢ = O, and the set 


att = —e given by Eq. (41), or simply the negative of 
the displacements that were found at ¢ = ¢«. In the 
actual case, no displacements are present at = —e, 


but these displacements may be regarded as being of a 
fictitious nature, the only purpose of which is to allow 
the step-by-step evaluation of the displacements to be 
started easily. With these sets of displacements, the 
next set may be obtained by application of Eq. (39). 
With the newly found set and the preceding two sets of 
displacements, the next set may then be found, etc., 
until a sufficient time history of the displacements is 
found from which maximum loading conditions may be 
determined. 

The reason for using the difference form of the deriva- 
tives as given by Eqs. (4) and (5) might now be given. 
Eq. (35) may be likened to the differential equation 
that relates the load to the deflection for a beam. The 
unknowns are the deflections at time m. The right- 
hand terms correspond to the loading, the first term 
being the only one that is not known since it contains 
the unknown deflection. The subsequent inversion of 
the matrix [D] leads to, in effect, the solution to this 
differential equation and, in the beam analogy, cor- 
responds to the integration of the loading to obtain the 
deflection. When numerical methods are used, the 
deflection may be computed with good accuracy by in- 
tegration of the loading. On the other hand, if the 
difference equations that apply at an interior ordinate 
had been used, the matrix [A] would have appeared as 
an operator on one of the known deflections on the 
right-hand side of the equation. Effectively, its opera- 
to differentiate a known deflection, 


tion would be 
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corresponding in the beam analogy to the attempt 
obtain the load that caused a given deflection. This 
process, however, cannot be done with accuracy whey 
numerical methods are used because of the difficyly 
encountered in the form of small differences of large 


numbers. The difference equations that apply at ay 


outer ordinate and, consequently, lead to an integratigy 


process therefore have to be used. 


COMPUTATION OF LOADS AND STRESSES 


Once the time history of the displacements has bee 
found from Eq. (39), the normal loading on the wing ca 
If the notation 
of Eq. (37) is used, Eq. (33) may be written 


p > [no] \w n + iO n (43 


This equation shows that the loading matrix /p, may bx 


be found with little additional work. 


found by adding an easily computed matrix to the 
matrix |Q!, the value of which will have already bee 
determined in the response calculation. 

The loadings thus found are considered to be applied 
statically, and the stresses are then found by ordinary 
static means. Since the loadings can be computed for 
any value of time, the complete stress history of any 


point in the structure may be computed. In general, 








the maximum stress at different points in the structure | 


is expected to occur at different times. Some guide 


as to the probable time of occurrence of the most severe | 


stress may be had, however, if the computed wing de- | 


flection is observed. It is likely that maximum stress 
occurs in the range where wing bending appears to be 
the most pronounced. 

The acceleration of any point in the structure may be 
found, if desired, with the aid of Eq. (5). 


EXAMPLE 


As an illustration of the method of analysis givet 
herein, the response of a two-engined commercial ait 
plane due to a sharp-edge gust is presented. The mass 
variation over the wing semispan, the station locations 
and the equivalent mass concentrations are shown ti 
Fig. 7. A solution is made for a forward velocity 
flight of about 210 m.p.h. and a gust velocity of 10 it 
per sec. Only the results of the solution are presented 

The time histories of the loads that occur at each | 
the spanwise stations are shown i1i Fig. 8. These curves, 
when judged on the basis that the half period of tht 
fundamental covers about 15 time intervals, seem t 
indicate the presence of some second-mode excitatiol 
in the response. The stresses that occur at stations! 


1, and 2 are shown in Fig. 9. 


CONCLUDING DISCUSSION 


A method for computing the stresses and structural 


action of an airplane flying through a gust has bee! 
given. 


The method is based on aerodynamic sttip | 
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eory, but approximate overall corrections for com- 


th 
pre 


ssibility and three-dimensional effects can be made, 
as iS indicated by a suggested correction procedure. 
Although the analysis presented is given in terms of 
airplane vertical motion and wing bending, it can be 
readily extended to the more general case where pitch, 
wing torsion, fuselage bending, and tail forces are also 
included. The result would be a recurrence equation 
similar in form to Eq. (39) but with larger order ma- 
trices. The form, for example, of the 7-matrices that 
would occur in |Q) [see Eq. (37)] for the more general 
case is illustrated by the matrix given in Fig. 10, where 
All ele- 


It may be of interest to 


the numbers shown are chosen at random. 
ments not shown are zero. 
explain briefly the significance of the various terms that 
appear in the matrix. In order to facilitate the explana- 
tion, the matrix has been partitioned into several sub- 
matrices. The terms in the upper left-hand box are 
associated with wing bending and the airplane vertical 
motion (the case treated herein), whereas the terms in 
the lower right-hand box are associated with wing 
torsion and airplane pitching. The terms along the two 
subdiagonals serve to couple together the bending and 
twisting action. The terms in the first row and first 
column are associated with fuselage bending. 

In a good many cases that may be considered, wing 
twisting and fuselage bending may prove to be of neg- 
ligible importance. Some investigators have indicated! 
that, unless the forward speed of the airplane approaches 
the flutter or divergence speed of the wing, the tor- 
Like- 
wise, in cases in which the fuselage is rather stiff, the 


sional deformations do not have to be included. 


effect of fuselage flexibility on the response may be 
rather small. In such cases in which either or both of 
these flexibilities may be ignored, the analysis is, of 
course, simplified and shortened. In the present state 
of understanding of gust response analysis, enough in 
formation is not available to indicate definitely when 
and when not to include the various flexibilities of the 
aircraft structure, and analysis of the type indicated 
herein may provide a useful means to assess their im- 
portance. 

Both the symmetrical and antisymmetrical types of 
gust can be handled by the analysis given in this paper 
In general, the symmetrical gust is expected to produce 
the most severe stress condition, and therefore only the 
flexibility matrix that applies to the symmetrical case 
has been given. This matrix was derived by using the 
boundary conditions for the symmetrical deformation 
of a free-free beam. The case of an antisymmetrical 
gust can be treated by replacing this matrix by the one 
that applies for the antisymmetrical deformation of a 
lree-free beam. The case of a general unsymmetrical 
gust can be handled by first breaking the gust into two 
parts, a symmetrical part and an antisymmetrical part, 
and then treating each part independently. 

It might be of interest to compare briefly the matrix 
method given in this paper to a modal function solu- 


OF 
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tion. One of the chief disadvantages of the modal 
function solution is that the modes and frequencies of 
natural vibration of the structure have to be computed 
in advance. Then a large number of integrals that in- 
volve these modes have to be determined in order to set 
up the problem. In the present analysis, the physical 
properties of the airplane are used directly in the setting 
up of the problem. Further, in order to make the modal 
solution practical, the higher modes must be dropped 
and only the basic or fundamental modes can be used. 
Hence, the success of the analysis depends to a large 
degree on how well single modal functions, one mode 
each for bending and torsion, can represent the airplane 
distortion. 
distortions are found, for all practical purposes, as the 


In the analysis of the present paper, the 


correct values at a number of spanwise stations, at 
least to within the accuracy to which the aerodynamic 
and structural parameters are known. 

One of the important features of the method pre 
sented is that it is not restricted to the gust problem. 
The approach may be used successfully to treat other 
problems of a similar nature. The landing problem 
can be handled by simply replacing the distributed 
In the 
landing problem also, the problem is less diflicult to set 


gust force by the concentrated landing forces. 


up, since the aerodynamic terms do not ordinarily have 
to 
which aerodynamic forces are included may be solved 
The approach 


be included. The landing problem, however, in 
by the method of this paper if desired. 
may also be used to solve the problem of the release of 
heavy objects such as bombs. This problem could be 
considered the inverse of the gust problem, a load 1s re 
leased rather than encountered. Some maneuvering 
problems, such as the sudden deflection of the ailerons, 
and a number of their transient problems might also be 
treated by an approach similar to that given in this 


paper. 


Appendix A 


DEFINITION OF MATRIXES USED IN ANALYSIS 


For convenience in presentation, most of the ma 
trixes and matrix elements that are used in the analysis 
are defined in this appendix. 

Matrixes.—The various matrixes that are used in the 
analysis are defined as follows for the case in which the 
wing semispan is divided into six sections: 


w(Q) p(O) 
w(1) pi) 
w(2) >(2) 
1 = (3)? P| = oe 
w(4) p(A) 
w(d) pl) 


A| 





See Appendix B for derivation of 
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ni(0) f = Belv (Ag 


ni(1) 
ni(2) g = Bcle®ye~ - (A5 





n(5) J 





[D] = [A] — [nl] Appendix B 








a a DERIVATION OF MATRIX EQUATION FOR EQUILIBRIVy 
| 42) F(2) OF LOAD 

n* f(3)) F| = F(3) In accordance with the assumptions made in this inw 

f(4) F(4) paper, the wing semispan is considered to be divided 

f(5) F(5) into six sections, with a station point at the center of 

4 each section (see Fig. 3). The inertia force of the mass 

80) and the aerodynamic force that develops over each gee. 

g(1) ™ tion are in turn assumed to be concentrated at the 

lg] = ate) ' respective station points. The wing semispan is thys 
g(3) 2(4) effectively a beam bending under six concentrated loads The 
z | and, as such, will have a linearly varying moment be. IS Z€1 

i go) tween each station. The wing is further assumed to | 


Matrix Elements.—The matrix elements that appear have a linear 1/// variation between each station, with 
' 


in the matrixes defined in the previous section are pre- the correct or designated value of EJ at the station | Fron 


sented for convenience in terms of the following common points. - atin 
factors: With these assumptions, the //E/ variation be. : 
i ietnciadlll, 4 y tween each station may be represented by an analytical 
eke Mei ie, Bye ae (Al) function. Application of the principles of engineering 
} beam theory then allows the exact calculation (for the 
in which the last four are associated with the @ function assumptions stated) of the deflection of each spanwise 
[see Eq. (23)]. With these factors, the elements that station relative to the wing centerline, where for the 
must be computed at each spanwise station are as fol- case of symmetrical deformations the slope must be 
lows: zero. From these deflections the deflections relative to 
r - station 0 may readily be determined. Although the whic! 
ad = — me (1 — ©o)Bcl + (« + 5 ty) el details of this calculation are not presented, it can be 
eee a (A2) shown that the deflections relative to station 0 may be . 
dy = —(3/26) (1 — &) Bel given by the following matrix equation: Subs 
oo ~ te \@(1) an a2.0 0 0 | |M(O) 
no = —(2m/e?) + do (2) b? (2, Ax Qe; 0 O M(1) Mult 
a= (Rete /e2 mys _ ‘ | [ 
m = (5m/e ) + d, (A3) (3) = EI(0) 31 Q32 433 A34 O we) (BI [7] 
nm = —(4m/e*) + dp w@(4) (41 As C43 O44 45) |M(3) 
nz = (m/e?) + ds (5) 51 As2 G53 M54 55| |M(4) 
This 
flecti 
case 
where the matrix elements are defined by the equations body 
Qi = Not + 7A, ) to se 
da = o(Ar + Az) + 2d + ABs rice 
31 = No(Ar + Ao + As) + AL7Ad + Ai(Ae + Az) Bi datur 
Gs, = do(\r + Az + As + As) + MZAL + Ar(Ae + As + Ag)Bi : ts 
ons 


51 = Ao(Ar + Ao + As + As + As) + AI7A1 + Ar(Az + Az + Ag + As) Bi 


(B2 show. 


aa = AVC, + AAD, + ho*Ao 

G32 = APC, + Ai(A2 + Az)Di + Av?Ao + AcdgBe 

dy = PCr + Ar(A2 + As + Ag)Di + Av*Ae + Ao(As + Ag) Be 

dsz = A7Ci + Ai(Az + Az + As + As)Di + Av?Ae + Ao(Azs + Ag + As) Be 
| Eq. (B2) continued on facing page | 
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ae = de? C2 

33 = Ao*C2 + AedgDe2 + Az7Az 

aig = h2?C2 + Ao(A3 + 4) De + \3°A 3 + A3A4Bs 

sg = Ao*C2 + Ao(Az + Ay + As) De + Az"Az + Az(Ay + As) Bs 

a3, = A3°Cz > (B2, Cont.) 


au = A32C3 + A3AyD3 + Ny"Ag 

= N3?C3 + Az(Aq + As)D3 + Ag7Ag + AASB 
as = NC, 

a5 = MC, + AsAsD, + 357A; 


54 


in which 


1 J(0) 1 J(0) 1 J(0O) 1 (0) 
A, ona P + < 2 B; = > ° 
4J(4 — 1) 12 J (1) 3/(4 -— 1) 6 I(t) 

(B3) 
ae | I(0) 1 J(O) —e 1 J(0O) 1 [(0) 
~ 4271-1) 127@’ ~° 616-1) ° 3I@ 


The moment ./; does not appear in Eq. (B1) because the boundary condition is used that the moment at station 5 
iszero. For convenience, Eq. (B1) may be given the abbreviated form 


@, = [b?/EI(O)] (i) | (B4) 


From the five loads p(1), p(2), p(3), p(4), and p(5), the moment at each station may be found. The equations re- 
lating the moments to the loads can be shown to be given by the matrix equation 


M(0) Mi Ar tAe Ar tAetAs Ar tA2+As+ AG Ar + Aet+ As + Ag + As ip(1) 
M(1) 0 Ne Ae + Az Ao + As + Ag No + Ag + Ag + As P(2)| 
M(2)' = 6| 0 0 A3 As + Aq As + Ag + As | /P(3) (B5) 
M(3) 0 0 0 Nj As + As | |P(4) 
M(4) 0 0 0 0 As_| |p(5) 
which can be abbreviated simply To aid in the derivation, [E/(0)/b*] [ [4] [He] ]~ is now 


written in the general form 


M\ = b [fi] |p (B6) 

bir Dia bis dis bis| 

. Dia bee beg bes bes| 

@, = [b*/EI(0)] [2h] [2] |p| (BY?) ” [ [21] [F2]]~* = [brs bes Bss bss bss] (B11) 
p 3 


Dis boy bay bay 45 


Substitution of this equation into Eq. (B4) gives 


Multiplication through by the reciprocal of [b*///(0)| i. a ie 
(| [H2] results in the equation 15 92% 935 945 65) 
Thus, with this equation, Eq. (B10) may be trans- 


formed in the form 


[E1(0)/b5} [ [2h] [2] ]~' |w| = |p| (BS) 


This equation thus gives the loads in terms of the de- w(0) 
® - . . . ° ee 

flection of each station relative to station 0. In the 

‘ . . . . boy by bie bis bis bis w(1) P( 1) 
case under consideration, however, the wing is a free a ee w(2) p(2) 

: : ° 02 2 Y22 U2 24 25 WZ o 4 
body capable of motion through space, and therefore, i ' “ae ‘ b ; b (3) = »(3) (B12) 

: “ . 203 O13 O23 033 D3 35 We ‘ 
to set up properly the equations of motion, the deflec- i rs b ‘i i b (4) p(4) 
. iis ; . % pci bos Dig bog 0; 1 ba5 | \w(- 
tion must be given relative to a fixed datum line. This eek ae a sae a 2 
bos bis 5 35 45 Bs5_} w(5) Plo) 


datum line is most conveniently located as the position 





of the wing prior to action of the disturbing loads. 
Consideration of the deflections shown in Fig. 11 will 
show, therefore, that the following relation must exist: 


b= w — w(0) (B9) 


2 


Substitution of this equation into Eq. (B8) gives 





EI(0) ' ; 
ps [ (2h) [2] ]~* jw — w(0)| = |p| (B10) Figure ||.- Notation for wing deflection. 











550 JOURNAL OF THE AERONAUTICAL 
where 

bu = —(bu + dis + dis + dis + dis) 

bor = —(di2 + bo + bo3 + be4 + bos) 

bos = — (diz + be3 + b33 + b34 + 35) (B13) 

bu = — (diy + bey + b34 + b44 + bas) 

bos = —(dis + bo5 + 535 + b45 + 055) 


Eq. (B12) is noted to express all the loads except p(0) 
in terms of the six deflections. An additional equation 
allowing p(0) to be expressed also in terms of the de- 
flections may be established by use of the condition that 
all the loads acting on the wing semispan must add to 


zero—that is, 


p(0) + p(1) + p(2) + p(3) + p(d) + p(s) = 0 
(B14) 


This condition automatically satisfies the two boundary 
conditions that the shear must be zero at the tip and 
centerline of the wing. Thus, if the five equations 
represented by Eq. (B12) are added and use is made of 
Eqs. (B13) and (B14), there results the following equa- 


tion: 


boow(0) + bo, w( i) + bogw(2) + bygw(3) + 
boyw(4) + bosw(5) = p(O) (B15) 
where 
boo = — (bo + boo + bos + bos + bos) (B16) 


This equation may now be combined with Eq. (C17) to 


give, finally, 





boo bor bor bog bos bos | 'w(O) p(O) 
bo bi bie bis bis bis | Ww) pl) 
bor bie bee box bey be5 | w(2) p(2) m= 
bos 61; b23 b33 b34 bes | w(3) ‘a p(3) wer) 
bos bis boy bas bay a5 | w(4) p(4) 
| bos O15 b25 35 45 Os5_] w(5) p(d) 





This equation is thus the desired matrix equation that 
relates the normal loads to the deflection. If the square 
matrix is denoted by [A], the equation may be abbre- 


viated conveniently to the form 
[A] w) = |p| (B18) 


which is the form used in the test, Eq. (13). 


Appendix C 


EVALUATION OF 


EXPONENTIAL 


EQUATION FOR THE 
INVOLVING 
KERNEL 


RECURRENCE 


DUHAMEL’'S INTEGRAL AN 


Consider the integral 


*/ 
aa_f wb(t — r)dr 


in which ® is given by Eq. (23). 


(23) into this integral gives 


(C1) 


Substitution of Eq. 


SCIENCES 


ae 
Baliye" f we dr 
0 


Mathematically, the integral in this expression m 
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ay 

interpreted to represent the area under the pate 
given as a product of w and e”’. In accordance with 
numerical evaluation processes, the interval 0) to ; maj 
be divided into a number of time stations of interval , 
The product of w and e”’ may then be found at each oj 
the time stations, and from these products the pe 
under the curve may be determined in first approxim, 
tion by the trapezoidal method of determining areas 


Thus, if the mth time station corresponds to time ¢, the 


value /, of expression (C2) may be approximated a 
| 


follows: 


i oe = Bcl@yee~ *"* [wy e"* of weer + ois 


n—1) ‘ yn 
wpe" "* + (1/2) w,e™] 


where wy does not appear, since the initial condition js 
used that w is zero at = 0. More accurate methods 
such as Simpson's method, could be used for deter 
mining the area under the curve, but, because of th 
small interval chosen, the consequent increase in ac 


curacy is negligible. If the notation 


F,, = Bcldyee~%"* [wie + wre” +... w,e™"~ 4] 


is introduced, Eq. (C3) may be written simply 


I, = F, + (1/2) Bclebow, C5 


Suppose now the expression (C2) is expanded similarly, 
only for an upper limit of ¢ — e, the expanded result 
would be 


¥ oN 1 5 72 
T,-1 = Bcl®yee~*'"~"* [wie™ + we™ +... 


Wn “ o— + ( ] z Wy, ie i " Ch) 
By analogy with Eq. (C4), however, 
F > Bcl®yee z : [we + wre’ ee 
w,2e'"~*} (Ci 


A study of Eqs. (C+) and (C7) will show that th 


following useful relation must therefore exist: 


F, = e “F,1 + Bcledye~ 2 (CS 


Wy —1 


Thus, Eq. (C5) when used in conjunction with this 
equation allows for the convenient step-by-step evalu 
(C1). Rather than having t 
carry through a complete evaluation for each successivt 


ation of expression 


value of ¢, the value of the expression at one time in 
terval can be obtained directly from the already estab 
lished value at the previous time interval. 
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On the Thickness of a Steady Shock Wave 


R. von MISES* 


Harvard t 


SUMMARY 


A simple formula is derived, not requiring any elaborate integra- 
tion, which gives exact and narrow limits for the extension of the 
transition layer in the one-dimensional, viscous, heat-conducting 


steady flow. The method employed also lends itself to a quick 


survey of all possible flow patterns. 


(I) THE Basic EQUATIONS 


Fu, p, p, AND 7 DENOTE the velocity, pressure, density, 
I and (absolute) temperature of a steady flow in x- 
direction and if m, C,, and C, are three constants, the 
following three equations hold: 


pu m 


mu+p—a Cm (2) 


gR 


m (“ - r) + u(p — o) —- H = Cam (3) 
2 ¥ ] 

Eq. (1) is the continuity condition; Eq. (2) is the mo- 
mentum equation, o being the tensile stress due to vis- 
cosity; Eq. (3) expresses the specifying condition that 
the total heat input for an element is given by the in- 
tensity // of the heat flow. It is here assumed that the 
fluid is a perfect gas—i.e., it fulfills the equation of 
state 


pb = gRpT (R = const.) (4) 
and the specific heat at constant volume is given by the 
factor of Tin Eq. (3). 

The usual assumptions for viscosity and heat flow 


are 


a = (4/3)u(du/dx), H = k(dT/dx) (5) 
The viscosity coefficient » and the conductivity k are 
known as functions of 7. Both increase with T at 
about the same rate, so that the quotient u/k remains 
sensibly constant. 

If Eqs. (4) and (5) are introduced in Egs. (1) to 
(3), the following first-order differential equations for 
wand 7 result: 


4u du IRT : 
. =ut e — Cj (6) 
3m dx u 
k dT : ; u? = ... a 
= ¢ — C, _ + 1 (4) 
m dx 2 eel 
Received March 2, 1950. 
* Gordon McKay Professor of Aerodynamics and Applied 
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(II) DIMENSIONLESS VARIABLES 
It is seen that C, has the dimension of a velocity and 
C, that of a velocity squared. We may therefore use, 
instead of u and 7, the dimensionless variables v and 
94 

vy = 4?/2C,? 
gRT p a* (8) 

C2 CYp) Cy? 


where a? = yp/p is the so-called sound velocity squared. 


It is also seen that 


2 = (2 


y) (v/t) (9) 


Introducing Eq. (8) in Eqs. (6) and (7), the equations 
become 


4 lv ‘ 
ee t — (V2v — 2v) (10) 
3 m dx 
y-—lkdt y 
={-— —l)(e—-vQ+e) (11) 
gR mdx ty \ — 
where c = C,/C,*. To eliminate x, we divide Eq. (11) 


by Eq. (10) and use the notation (Prandtl Number) 


P = [y/(v — 1)]gR(u/k) (12) 
This leads to 
dt _ 4P t —_ to (13) 
dv 3yt—t, 
with /) and ¢,, two given functions of v, 
to = (y-— 1) @ — V+ o)\ 
/ ; (14) 
t., = V20 — 2D { 


The first-order differential equation [Eqs. (13) and 
(14)] has for given c a simple infinity of solutions. 
Once the relation between ¢ and v is found, Eq. (10) 
gives v as function of x and Eq. (11) gives ¢ as function of 
Except for scale factors and 
different flow patterns, 


x by simple quadratures. 
the origin of x, there are ~? 
«© ' for each c-value. 


) 


(III) THe SINGULAR Points OF Eg. (13) 


Whether the Prandtl Number is constant or not, if 
only 


0< P<-« (15) 


Eq. (13) has a singularity wherever the numerator and 
denominator of Eq. (13) vanish simultaneously. The 


conditions / = ¢,, and¢ = fy lead to 
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ty / 
t —_— JO = D9), 
\ ‘20, aa ~~" | 
¥ (19) 
. 4+ uy = te + v2 | 
y-1 y¥- 1 











FIG. | 


(16) 


The first of these equations represents in the v,/-plane 
a parabola, not dependent on c, tangential to the ¢-axis 
inv = ¢ = O and passing throughv = '/2, t = 0 (Fig. 1). 
The second Eq. (16) is that of a straight line of constant 
negative slope, its location depending on c. The 
two points of intersection, 1 and 2, are the singular 
points. 
Elimination of ¢ in the Eqs. (16) gives 


(y + l)v -— yV 2Qv + (vy —l)ce =0 (17) 


This equation has two distinct roots, 2, v2, both corre- 
sponding to positive ¢ if and only if 
49 
48 


l 7? 
<< = 
2 2(y? — 1) 


(18) 


At the upper limit 7; and v7, coincide; at the lower limit 
one v takes the value '/. and t becomes zero (dotted 
lines in Fig. 1). In what follows we shall assume that 
Eq. (18) is fulfilled. 

If the greater of the two roots of Eq. (17) is called 
2%, the well-known methods of geometrical analysis 
show, for Eq. (13), that point 1 is a nodal point and 
point 2 is a saddle point. Exactly two integral curves 
pass through point 2 and an infinite number through 
point 1, all of them except one with the same tangent. 
There is one and only one integral curve passing 
through both singular points. This integral corre- 
sponds to the transition between the two states repre- 
sented by 7, f; and vs, f. Since in the region between 
the curves ¢ = ¢,, and t = fy (also indicated in Fig. 1) 
t—t. << Oandt — ty > 0, it follows from Eqs. (10) and 
(11) that, for the transition line, dv/dx is negative and 
dt/dx is positive. In other words, the transition goes 
from larger v (smaller p) and smaller ¢ toward smaller 
v (larger p) and larger /. 

On the other hand, it is seen from the two Eqs. (16) 
(both of which are fulfilled by 7, 4 as well as by 72, f2) 
that 


If these are rewritten in the original variables 
and ¢ by p/ pC," 


that is, 


if v is replaced by u?/2C,? they read 


pi + muy, = po + Muy 
Y Pr uy" Y p2 , Us" 
%¥~- in 2 y—Il1lp 2 


Thus we have found, under the sole condition (15), that 
the states corresponding to the singular points | and? 
fulfill the Rankine-Hugoniot conditions and that the 
transition always implies a compression. 

According to Eq. (9) the lines of constant Mach 
Number in the 2v,f-plane are the rays through the 
origin. In particular, the line J/° = 1, orv ft = y,2, 
passes through the above-mentioned point on the para- 
bola t = ¢,. where the two singularities coincide. It 
follows that each transition starts in a supersonic and 
ends in a subsonic region. 

Since 7%, ¢; and v, tf, fulfill Eqs. (16), it follows from 


Eq. (9) and the first Eq. (16) that 
I 4h 1 ; 
Y M;? 20, V/ 20; | 20) 
ts ' “= 
Y AM, 7 2v V/ 20» 


One can therefore express the four coordinates in terms 
of M, and M,. If these expressions are introduced in 
the second Eq. (16) and if c is eliminated, the well- 
known relation 


2yMi2M2? = 2 + (y — 1) (Mi? + M2) 21) 


results. 


(IV) THE INTEGRAL CURVES IN THE ?,f-PLANE 


For any constant Prandtl Number P, it is easy to 
plot the integral curves in the v,f-plane by means of the 
isoclines, even without any previous knowledge of the 
nature of the singularities, etc. 

An isocline of the differential Eq. (1: 
P = const., is a curve on which the ratio 


3), in the case 


(t — to)/(t — t,) 


has a constant value. Let us call this ratio A, so that 
Eq. (13) reads 

dt/dv = 22 
and ¢, the ordinate of the curve corresponding to a col 


stant A 


ty — to to acy, A (9 
tr —t 1—A 
It is seen that ¢/,, corresponds to \ = © and fo corre 


(4P/3y)d (22) | 
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Integrals for M,=2 








sponds to = 0. All isoclines except those for \ = 1 
and \ = —(y — 1) are ordinary parabolas passing 
through points | and 2 and tangential to the f-axis. For 
\ = 1, the parabola degenerates into the parallels 
the 


to the axis through points 1 and 2. In case 


\ = —(y — 1), the second Eq. (20) gives 


th, = (c — v) (y — 1)/y¥ 


that is, the straight line connecting points | and 2. 

Fig. 2 shows the set of isoclines for 14, = 2, valid for 
all constant values of P, and a set of integral curves 
for P = 1. For each of the latter curves giving ¢ as 
function of v or v as function of ¢, one can find the rela- 
tion between x and v and that between ¢ and v by 
quadratures from Eqs. (10) and (11), 


kdt 


| *v  ydy y-1 “t 
3m t — to gRm t 
The differences ¢ — tf and t — ¢ 
the drawing and allow an immediate qualitative dis- 


s = (24) 


appear immediately in 
cussion. 

(V) THE TRANSITION CURVE 
The main interest belongs to the singular integral 


passing through points 1 and 2. dt /dv, 
under which an integral passes through one of the 


The slope t’ = 


singular points, is found by applying the L’Hospital 
tule to the expression 0/0 on the right-hand side of Eq. 
(13). One finds from Eqs. (13) and (14) 

. $Pt’ — (y — 1) [1 — C/V 22)] : 
i= (25) 


3 ho + 2 — (1/v/ 20) 


Fory = v and for v = v», this is a quadratic equation 


lor’. It has, in each case, two real roots—one posi- 


A STEADY 


SHOCK WAVE 


tive and one negative. The negative ones—they may 


be called ¢,;’ and t,’ transition 


(20), Eq. 


are the slopes of the 
curve in the two end points. Using Eqs. 


(25) can also be written 


‘ , l 4P iy — 1)P 
ete’. — — — “ — = Q (26) 
7M? 3Y 3y7°M 
where MJ = M, fort,’ and M = Mefor to’. 
In the special case P = * 4, these equations are ful- 
filled by 
ty’ = ty’ = —[(y — 1) | (P = 4) (27) 
independently of \/. This indicates—-what could 
have been pointed out much earlier—that for P = 


‘/; the transition curve is the straight line between 
points | and 2. In terms of the original variables u, 7, 
the equation of this straight line [second Eq. (16)] is 

u* Y = 

= fe gRT = 

2 y¥-I!1 


const. (P = 3/4) 28) 
a well-known result first found by Becker. ' 

In the general case, Eq. (25) expresses the fact that 
the slope under which the transition curve begins in 
point | or ends in point 2 equals the slope of the isocline 
on which ¢’ has the value ¢,’ or f2’, respectively, and for 
which d is, according to Eq. (22), 


(3¥ 4P)t,’ or (3 4P)t’ (29) 


On the other hand, it is seen that the transition line 
cannot cross either one of these two isoclines. Thus, if 
a is the greater and 6 the smaller of the two (negative) 
values (29), one has the decisive relation 

on oe (30) 


a f 


In the case P = */4, the upper and lower limits coin- 


cide. For P + 3/4, itis seen that 

SY ; 3Y ca Aes 3 

a= t.", B= t f# P< 
4P 1p ' | 
(31) 
3Y SY os 3 

eo t “ B= to’ if P > 

{Pp ' 1P 

(VI) Tuer TRANSITION WITH RESPECT TO 4 


If Eq. (30) is introduced in the differential Eq. (10) 
and if Eqs. (14) and (23) are used, it follows that 
4yu dv lb —t 


ei wes < tti~t. = 
3m dx l1—BpB 


Since yu is increasing with temperature, one has 


Mm OuS be (33) 
where yw; and pw are the values of uw at the end points 
Therefore, from Eq. (32), 
3m dx dv 3m dx ? 
Ee < Sie (34) 


4, | - @ a = — Ie dy | — 6 
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ee 


(y + 1) (We -— Vin) (Wo — Vn) 


In the integral formula 


/ dv 
7 ] / ; 
(Vv — Vu) (V0 — V2) 


each ~/v can be replaced by u, since the scale factor 
C,v/2 drops out. If the integral is taken in the limits 
u’, u" with 


u’ = um — €(u, — Ue) | ae 
4 < l 2) 3¢ ) 
u" = Us + €(u, — Ue) f te ; - 
its value is 
l \/ur + ues a 
2 tog ( - 1)(* )ar (37) 
€ Uy — Us 


and the inequalities [Eq. (34)] give, for the correspond- 


ing interval L = x” — x’, 


4ml—a 


D < tml —B 
3smyt+ 1 


<7 & 
3myt+ 1 


I (38) 


Let us call Ly the limiting value that L assumes for 
some constant u = yu* when heat conduction is entirely 
neglected—that is, P = ~,a = B = 0. Then it fol- 
lows from Eqs. (37) and (38) that 


2 4 l * fa U2 
Lo = log ( _ 1) i ( + ) (39) 
¥+13 € m \t — Ue 


and instead of Eqs. (37) and (38) one can write 


tt «we, £4 = eo OL, (40) 


u* mM 


* 

All one has to do, once Ly is known, in order to find 
the limits of L is to compute the negative roots /,’ and 
t,’ of the quadratic Eq. (26) (the first for 1J = As, the 
second for \J = M,.); to determine a, 8 from Eqs. (31); 
and then to insert these values in Eq. (40). 


(VII) EXAMPLE 
Let us assume the initial Mach Number J/; = 2. 

From Eg. (21) with y = 1.4 follows MJ.” = '/3. The 
negative roots of Eq. (26) are 

t,’ = —(. 13, t,’ = —().22, for P = 4 2 

ty’ = —().20, to’ = — 0.34, ar = |] 
and from Eqs. (31) 

L—e@= 1.46, 1— 8 = 1:88, for? = '/, 

l—ea io. t= B82 tas foe r= ! 


To find an estimate for Lo, we call p* the density at 
that point of the flow where 


u = (% + Ue)/2 


2 [ve log (W011 — V0) — Vive log (Wu — vo) 
=2 om (35 


‘a / 
VU — Vit2 


Then (#; + ue) in Eq. (39) can be replaced by 2m, 
Taking « = 0.05, the length over which 90 per cent oj 
the velocity drop occurs is 


Lo = 5.6 (u*/p*) [1/(u — u2)] 


in good agreement with the approximate result of Tay 
lor.2. With an average value for u*/p* ~ 1.6 X 10: 
sq.ft. per sec. and u; — wu, = 10 ft. per sec., the lengt, 
Lo will be 0.9 X 10-4 ft. = 0.03 mm. This gives, j 
the variation of yu is neglected, 


for P = '/, 
for P = ] 


0.044 < L < 0.056 m/m, 
0.036 < L < 0.040 m/m, 





In Fig. 3 the situation in the v,/-plane for M, = 2s 
shown. The ordinates are magnified at the rate 5:] 
The solid lines are the isoclines; the broken lines are the | 
integral curves for P = '/o, */4, and 1. 

Recently, the case of P = 0.75 has been studied, with 
essentially the same results, in the paper by Morduchow 
and Libby.* In an extensive study by Reissner and 
Meyerhoff,! some special types of motion were analyzed 
by numerical methods of integration, starting from the 
second-order differential equation, which can be de- 
rived from Eqs. (6) and (7) by eliminating either T or 


u. 
(Continued on page 594) 
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The transition curves for P=0.5, 0.75, !. (M,* 2) 
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A Type of Lifting Rotor with Inherent 
Stability 
KURT HOHENEMSER* 
McDonnell Atreraft Corporation 


SUMMARY 


The more important of the forward flight aerodynamic char- 
acteristics of a rotor type are presented where the collective blade 
pitch angle is kinematically coupled with the collective flapping 
or coning angle, while no coupling is provided between cyclic 
flapping and cyclic pitch angles. The rotor with pitch-cone 
change may be designed to have positive static stability as com- 
pared to the negative static stability of the conventional fixed- 
pitch rotor; partial blade stall during pull-ups or upward gusts 
may be avoided in spite of a high and economic thrust coefficient 
in steady flight, and no adjustment of the collective pitch con 
trol is required for the transition from powered flight to autorota- 
tion. The rotor with pitch-cone change, because of its inherent 
stability, promises, without recourse to external stabilizing de- 
vices, appreciably improved flying qualities of the helicopter. 


DISCUSSION 


[' IS GENERALLY KNOWN that the conventional lifting 
rotor with hinged or tetering blades has character- 
istics that produce, in forward flight of the helicopter, 
longitudinal instability unless external stabilizing sur- 
faces or other stabilizing devices are provided.! * ® 7 
We first discuss the main parameters responsible for the 
longitudinal control and stability characteristics of the 
helicopter. 

Fig. | shows a schematic side view of a lifting rotor 
with hinged or tetering blades. JV is the direction of 
the relative air speed; a measures the attitude of the 
control plane or plane of zero cyclic feathering with 
respect to the direction of flow, positive if inclined 
backwards; and a, measures the backward inclination 
of the tip path plane with respect to the control plane. 
The thrust vector 7 is approximately perpendicular to 
the tip path plane and, in the undisturbed flight condi- 
tion, is assumed to pass through the aircraft center of 
gravity. x and z are body fixed reference axes through 
the ¢.g., « pointing into the direction of undisturbed 
flight and z perpendicular to x downwards. Changes of 
forward flight speed and of attitude produce the follow- 
ing three pitching moments about the c.g.: 

An increase in forward flight speed V’ without change 
of attitude a@ causes an increase in a, or a backward 
inclination of the thrust vector 7 producing a_ tail- 
heavy pitching moment about the c.g. __In other words, 
the lifting rotor has positive velocity stability. 

A slow increase in attitude angle a without change of 
flight speed |’ also causes an increase in d@; or a back- 

Presented at the Rotating Wing Aircraft Session, Eighteenth 
Annual Meeting, I.A.S., New York, January 23-26, 1950 

* Chief of Aerodynamics, Helicopter Division 


ward inclination of the thrust vector 7° relative to the 
fuselage producing a tail-heavy pitching moment about 
the c.g. In other words, the lifting rotor is statically 
unstable. The fact that the lifting rotor is unstable 
with change of attitude was obscured for some time by 
the stable stick travel observed in most helicopter types, 
which means that an increase in steady flight speed re- 
quires a more forward position of the control stick. 
For the airplane, a stable stick travel is an indication of 
static stability. This is not true for the helicopter, 
where a stable stick travel may be produced by the 
velocity stability in spite of instability with change of 
attitude. 

A rapid increase in attitude angle a without change 
of flight speed I’ causes a decrease in a, or a forward 
inclination of the thrust vector 7’ relative to the fuse- 
lage, producing a nose-heavy pitching moment about 
the c.g. proportional to the rate of change of attitude. 
In other words, the /ifting rotor has positive pitch damping. 

Finally, an increase in attitude angle a produces an 
The increase in attitude angle a 
downward 


increase in thrust 7. 


may also be described as an additional 
velocity component of the helicopter 2, and we may 
state that the lifting rotor has positive vertical damping. 
Although there is a total of nine longitudinal stability 
derivatives, numerical the complete 


equations of motion has shown that only the four de- 


evaluation of 


rivatives just mentioned are of major influence on the 
flying qualities of the helicopter. The effect of these 
four derivatives is the following: 

The velocity stability is, as was mentioned before, 
mainly responsible for a stable stick travel. The 
velocity stability should, however, be kept in moderate 
limits, because the frequency of the long-period oscilla- 
tion of the helicopter increases with increased velocity 
stability, and a low frequency is desirable from a point 


of view of flying characteristics. 
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Schematic side view of lifting rotor 
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Plan view of rotor hub with pitch-cone change. 


Fic. 2. 


Static stability is necessary for damping the long- 
period oscillation. The instability of most helicopters 
is not obvious under normal flight conditions, because 
the period of the instable oscillation is so long that the 
pilot has ample time to correct any tendencies of self- 
excitation. If, however, because of a wrong c.g. loca- 
tion or a lack of forward stick travel, the pilot encoun- 
ters a condition where he runs out of forward longi- 
tudinal control, the instability becomes obvious, and 
dangerous flight conditions with extreme and uncontrol- 
lable attitude changes occur. It may be noted here 
that a pitch governor keeping the rotor speed constant 
increases the static instability of the helicopter, and for 
this reason such a governor should be avoided. 

The pitch damping, which is proportional to the blade 
moment of inertia, has a marked effect on the flying 
characteristics of the helicopter. The higher the pitch 
damping, the longer and more damped is the oscillation 
period of the helicopter. Therefore, lifting rotors with 
either heavy blades or heavy tip weights of jet engines 
or rotor types, with artificially increased pitch damping 
like the Bell and Hiller rotors, require less attention by 
the pilot than conventional rotor types. An increase 
in pitch damping, however, slows down the control re- 
sponse, and the pitch damping therefore should not be 
excessive. 

The vertical damping is responsible for the load fac- 
tors produced by gusts. A low vertical damping is de- 
sirable from a point of view of smooth flight and also 
of good dynamic stability. 

One method of correcting the static instability of 
the conventional lifting rotor is to compensate it 
by the static stability produced 
stabilizing surface or by appropriate 
shape. In spite of difficulties caused by the rotor 
downwash, _ satisfactory been 
achieved in some cases by stabilizing surfaces at least 
However, the better method 


positive by a 
an fuselage 


stabilization has 
for a limited speed range. 
seems to be the elimination of the unstable moments 
at the source rather than their compensation by an addi- 
tional device. 

Another method of correcting the static instability 
of the conventional lifting rotor would be the use of 
large 6; angles for the horizontal blade hinge axes. By 
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this method, a, could be almost completely eliminated 
and the tip path plane would always coincide with the 
control plane. This method would, however, not only 
eliminate the static instability but also eliminate the 
two other pitching moment derivatives—the velocity 
stability and the pitch damping—both of which ar 
extremely important, as has been indicated before, 

A third method, which is the main subject of this 
paper, not only leaves the velocity stability and the 
pitch damping unimpaired but also provides positive 
static stability for the lifting rotor. The method cop. 
sists in coupling the collective blade pitch angle with 
the collective blade flapping or coning angle, while no 
coupling is provided between cyclic pitch and cyelic 
flapping angles. One manner of realizing this method js 
shown in Fig. 2, which is a plane view of a rotor hub 
with pitch-cone change. The hub is_ universally 
hinged at the rotor shaft, and the blades are individu. 
ally hinged to the hub at a distance from the rotor cen. 
ter. The blade pitch arm extends opposite to the direc. 
tion of rotation to the rotor center and is operated ina 
conventional manner by vertical links attached to the 
points A. 

The tip path plane of this rotor may be inclined by 
an angle a, with respect to the control plane without 
affecting the blade pitch. Any increase in coning angle 
do, however, produces a decrease in blade pitch angle 
proportional to tan 63. The complication of the hub 
and hinge assembly by having a universal center hinge 
and additional outer flap hinges is, as compared to the 
conventional rotor, more than offset by the elimination of 
drag hinges and their dampers, which are unnecessary 
in this design, and by the compactness and simplicity 
of the control mechanism. In case of a jet-driven 
rotor, the center universal hinge may be replaced by a 
self-aligning bearing, and the design, shown in Fig. 2, 
lends itself especially well to jet-driven lifting rotors. 

The following comparison between the conventional 
fixed-pitch rotor and the rotor with pitch cone-change 
is made on the basis of a simplified analysis following 
Bailey's N.A.C.A. Report® but using a constant profile 
drag coefficient of the blade of 0.012. The rotor 1s 
assumed to have rectangular untwisted blades with a 





~_— 





blade solidity of 0.06 and a blade inertia coefficient ol | 


10. 

Fig. 3, upper portion, shows, for an advance ratio y 0! 
0.35, the coning angle ap against attitude angle a of the 
control plane for different blade pitch angles 6. 

Fig. 3, lower portion, shows the backward tip pati 
plane inclination a; against a. For any given blade 
pitch angle 6, the backward tip path plane inclination 
increases with a, causing static instability of the lifting 
rotor. 

Now assume, as an example, that blade pitch angle 


and coning angle dy are related by the equation 
@ = 24 


—_ yee ao 


which is valid for a 6; angle of 70° and for a coming 
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Fic. 3. Coning angle a) and tip path plane inclination a, 


against a with blade pitch angle @ as parameter for advance ratio 
u = 0.35. Dashed lines for sample rotor with pitch-cone change. 





angle dy) of 5.5° belonging to a pitch angle @ of 9°. 
Points corresponding to this equation are connected in 
Fig. 3 by dashed lines. The higher the coning angle 
a, the lower the pitch angle 6. The dashed line in the 
lower part of Fig. 3 gives the change of a, with a, and 
it is seen that the slope is reversed as compared to the 
conventional rotor with fixed blade pitch angle 6. The 
rotor with pitch-cone change may be designed so that 
it provides positive static stability. 

Fig. 4, upper portion, shows the thrust coeflicient ¢, 
against a for different blade pitch angles @ and for the 
same advance ratio as before, w» = 0.35. The dashed 
line belongs to the sample rotor with pitch-cone change. 
It is seen that the increase of thrust with increased atti- 
tude angle a is considerably reduced as compared to the 
rotor with fixed blade pitch. 

Fig. 4, lower portion, shows the torque coefficient of a 
lifting rotor against a for different blade pitch angles 6. 
The dashed line, belonging to the rotor with pitch-cone 
change, does not show the curved character of the fixed- 
pitch lines and has a fairly constant slope. 

Diagrams similar to those presented here for an ad- 
vance ratio of uh = 0.35, were also developed for other 
advance ratios. The result is a complete comparison 
between the fixed-pitch rotor and the rotor with pitch- 
cone change, as shown in the following diagrams. The 
three quantities—backward inclination of tip path 
plane with respect to control plane a), thrust coefficient 
Cr, and torque coefficient cg—are plotted for both rotor 
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Fic. 4. Thrust coefficient cr and torque coefficient cg against 


a with blade pitch angle 6 as parameter for advance ration» = 


0.35. Dashed lines for sample rotor with pitch-cone change. 
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Fic. 5. Comparison between fixed-pitch rotor (above) and 
rotor with pitch-cone change (below). Tip path plane inclina 


tion a; against a with advance ratio uw as parameter 
lines for constant cg/cr equal to cg/cr in hovering. 
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types against a with » as parameter and against wu with 
a as parameter. The upper part of each figure con- 
tains the fixed-pitch rotor characteristics, assuming a 
pitch angle of @ = 9°; the lower part of each figure con- 


tains the characteristics of the rotor with pitch-cone 


change. 

Fig. 5 shows the comparison for the tip path plane 
inclination a; against a with w as parameter. This is 
the most important part of the comparison, since it 
shows the reversal of the slopes of the a; — a curves for 
ally values. The fixed-pitch rotor is statically instable 
in the whole range; the rotor with pitch-cone change is 
statically stable under all conditions. 

Fig. 6 shows the same characteristics in a different 
presentation namely, a; against » with a as parame- 
ter. The main character of the a; — yu curves is the 
same in both cases, indicating that the velocity stability 
of both rotor types will*be of the same order of magni- 
tude. 

Fig. 7 shows the comparison for the thrust coefficient 
cpagainst with was parameter. The slope of all cr — 
a curves is appreciably smaller for the rotor with pitch- 
cone change than for the fixed-pitch rotor. 

Fig. 8 shows the same characteristics in a different 
presentation—namely, cr against u with a as parameter. 
Fig. 9 shows the comparison for the torque coefficient 


(g against a with uw as parameter. 
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Fic. 10. Same as Fig. 9, except with uw as abscissa and a@ as 


parameter. 
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Fic. 11. Vertical damping 

against advance ratio u for constant cg/cr equal to ce/cr in hover 

ing. 1 for fixed pitch, @ = 9°; 2 for sample rotor with pitch 

cone change, @ = 24 — 2.7ao. Solid curves for constant torque; 
dashed curves for constant rotor speed 





— Z; and velocity stability Mj 


Fig. 10 gives the same set of values in a different 
presentation—namely, cg against u with a as parameter. 

The data given in these diagrams are sufficient for 
establishing a theoretical comparison between the flying 
characteristics of the fixed-pitch rotor and those of the 
rotor with pitch-cone change. The main improvement 
accomplished by the pitch-cone change is the positive 
static stability instead of the negative static stability 
of the fixed-pitch rotor. But there are also other im 
provements, the most important of which is the elimina- 
tion of partial blade stall at high forward speeds, es- 
pecially during maneuvers. Partial blade stall has 
proved to be a severe limitation of the lifting rotor, and 
it is accompanied by heavy vibrations and reduction or 
The main parameter controlling the 
Performancewise, 


loss of control. 
blade stall is the thrust coefficient. 
it is most economic to operate at a high thrust coeffi- 
cient and close to the condition of beginning blade 
stall.‘ In this case, however, pull-ups and upward gusts 
cause large portions of the blades to stall. The rotor 
with pitch-cone change reaches, during maneuvers and 
gusts, considerable lower thrust coefficients than the 
fixed-pitch rotor (see Fig. 7), and blade stall may there- 
fore be avoided in spite of a high and economic thrust 
coefficient in steady flight. 
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Fic. 12. Static stability — M; against advance ratio u for same 


conditions as in Fig. 11. 


A further improvement of the rotor with pitch-cone 
change as compared to the fixed-pitch rotor is its capa- 
bility of autorotation without adjustment of the collec- 
tive pitch control, which is a desirable safety feature. 

Finally, a quantitative comparison of stability deriva- 
tives of the two rotor types is of interest. Fig. 11 
shows the vertical damping Z; and the velocity sta- 
bility 7; against u for conditions with a power input 
equal to the hovering power. The solid lines have 
been computed under the assumption of constant torque 
and varying rotor speed, which is practically fulfilled in 
the case of slow oscillations of a helicopter without rotor 
speed governor. The dashed lines have been computed 
under the assumption of constant rotor speed and vary- 
ing torque, which is approximately fulfilled in the case 
of rapid motions of the helicopter when the inertia of 
the rotor prevents corresponding rapid changes of the 
rotor speed or in the case of a rotor speed governor con- 
trolling the engine output. Curves 1 belong to the 
fixed-pitch rotor, and curves 2 belong to the rotor with 
pitch-cone change. The vertical damping Z; is nearly 
the same in both cases if the torque is kept constant. 
This indicates that, for relatively slow pull-ups or in 
turns, the two rotor types will have the same vertical 
load factors. For constant rotor speed, however, the 
rotor with pitch-cone change has a considerably reduced 
vertical damping Z;, indicating a reduced sensitivity to 
sudden gusts. 

The velocity stability of the rotor with pitch-cone 
change is increased as compared to the fixed-pitch 
rotor, indicating a greater forward stick travel with in- 
creased forward speed. The increased velocity stabil- 
ity produces an increased frequency of the helicopter 
oscillation, which is not desirable. The 6; angle should, 
therefore, stay in the limits required for achieving static 
stability and should not be unnecessarily high. 

The pitch damping has not been plotted because it is 
the same for both types of rotor. 

Fig. 12 shows the static stability derivative J, for 
both rotor types, also against u. Positive M: is equiva- 
lent to static instability. The improvement by the 
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rotor with pitch-cone change is especially great in th 
case of constant torque (solid lines). 

Static stability is a necessary but not sufficient re. 
quirement for dynamic stability. Whether or not th 
helicopter is dynamically stable depends mainly on th. 
Fig. 13 indicates the 
maximum moment of inertia allowable in order to op. 


fuselage moment of inertia. 


tain dynamic stability with the rotor with pitch-cone 
change. The solid line represents the case of constant 
torque; the dashed line represents the case of constant 
r.p.m. In spite of the smaller static stability, the 
maximum allowable moment of inertia is larger in the 
case of constant r.p.m. than in the case of constant 
torque. If the inertia of the fuselage stays within thes 
limits, the helicopter with pitch-cone change will be 
longitudinally stable, without any stabilizing surface or 
stabilizing device, only by its own inherent stability, 
Fig. 13 indicates that the rotor stability increases with 
forward speed. Of course, the fuselage must have a 
shape so that its pitching moment derivatives do not 
offset the stability derivatives of the rotor proper, 
But even if the inherent stability of the rotor with 
pitch-cone change is not sufficient to overcome desta- 
bilizing effects of the fuselage, the use of this rotor type 
seems to promise appreciably improved flying qualities 
of the helicopter. 

To sum up the main characteristics of the rotor with 
pitch-cone change: 

(1) Positive static stability instead of negative static 
stability of the fixed-pitch rotor. 

(2) Possibility of avoiding partial blade stall during 
pull-ups or upward gusts in spite of a high and eco- 
nomic thrust coefficient in steady flight. 

(3) Elimination of collective pitch adjustment for 
the transition from powered flight to autorotation. 

(4) Considerable mitigation of loads caused by sud- 
den gusts. 

(5) Increase in velocity stability. 

(6) Unimpaired damping in pitch of the helicopter. 
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Fic. 13. Maximum allowable helicopter moment of inertia 
about pitching axis against » for dynamic stability when using 
sample rotor with pitch-cone change. Solid lines for constant 
torque; dashed lines for constant rotor speed. 
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(7) Unimpaired vertical maneuvering load factors if 


th 


of 


e rotor speed is allowed to vary without interference 


a rotor speed governor. 


(8) Possibility of attaining dynamic longitudinal sta- 
pility in forward flight without stabilizing surfaces or 


other stabilizing devices. 


We were mainly concerned here with the longitu- 
dinal stability of the rotor with pitch-cone change. 
Other problems related to this type of rotor (such as 


lateral st 


ibility, rotor speed oscillations, and coning 


angle oscillations) have been investigated and will be 


presented separately. 


SYMBOLS 


velocity of undisturbed flight 

attitude angle of control plane (plane of zero cyclic 
pitch) with.respect to flight velocity; positive if 
helicopter is inclined backwards 

angle of backward inclination of the tip path plane 
with respect to the control plane 

thrust vector, assumed to be perpendicular to the tip 
path plane 

body fixed reference axes, x pointing in the direction 
of undisturbed flight and z perpendicular to x 
downwards 

collective flapping or coning angle 

angle between outer kinematic flapping axis and tan 
gential direction 

profile drag coefficient of blade 

cpaR*/Iz = blade inertia coefficient 

blade chord 

air density 

slope of blade lift coefficient against blade angle of 
attack (assumed to be 5.7) 

rotor radius 

blade moment of inertia about horizontal hinge 

bc/rR = blade solidity 

number of blades 

blade pitch angle 

T/pQ?xR* = 

angular rotor speed 

Q/pQ?2xR® = 

rotor torque; positive if engine puts power into rotor 


thrust coefficient 
torque coefficient 


V cos a/Q“R = advance ratio 

speed of axial flow through rotor 

functions of » tabulated in Table 1 

small disturbance velocities in direction of x axis and 
of 2 axis 

small disturbance angular velocity about pitching 
axis 

gross weight of helicopter 

acceleration of gravity 

force in direction of z axis 

pitching moment; positive if tail heavy 

helicopter moment of inertia about pitching axis 

OZ /0z 

OM /ort 


(w/g)/pRA = mass coefficient 


(derivatives 


rotor disc area 

W/(V2/2)pA =~ 2Cr/p? = rotor lift coefficient 
time 

M;/ARpv = velocity stability 


—M3/ARpv = static stability 
—M,/AR?’pv = pitch damping 
—Z3/Apv = vertical damping 
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~ 


= I/(w/g) R* = dimensionless helicopter moment of 


inertia about pitching axis 


1/E = coefficients of frequency equations 
Ic = Q/AR(p/2)v? = 2cg/u? = torque coefficient 
qe = Og-/O : ‘ 
wen Ode/ CH \ derivatives 
dea = Og-/0a 
CM = M/AR(p/2)V? = pitching moment coefficient 
Cup = OC /Ou 
OC > derivatives 
CMa = 
Oa / 
h = distance between rotor hub and helicopter c.g. 


h = h/R= 


dimensionless distance 


AERODYNAMIC ROTOR CHARACTERISTICS 


According to reference 3, the following equations are 
available for determining the thrust coefficient cr, the 
torque coefficient cg, the backward inclination of the 
tip path plane with respect to the control plane a, and 
the coning angle a as functions of the three independent 
variables a, uw, and 6: 


2u*?), A= - V cr 2 for hovering) 
2cr/ao = A\ + BO 

2cg/ao = (C/a\b — Dd? — EXO — FO? 

a, = Gr + H0 

ao/y = IX + KO 


a = (A/p) + (er 


The rotor is assumed to have untwisted rectangular 
blades with a blade solidity of ¢ = 0.06, with a lift coeffi- 
cient slope of a = 5.7, with a constant profile drag coeffi- 
cient of 6 = 0.012, and with a blade inertia coefficient 
al y = 10. K as 
functions of uw are given in reference 3 for a tip-loss fac- 


The values for the coefficients A ... 


tor of 0.97 and are summarized in Table 1. 


TABLE | 
m 0 0.15 0.25 0.35 0.45 
A 0.470 0.476 0.487 0.504 0.529 
B 0.300 0.315 0.333 0.361 0.399 
: 0.250 0.256 0.266 0. 280 0.299 
D 0.470 0.502 0.558 0.646 0.766 
I 0.300 0.369 0.490 0.689 0.986 
F 0 0.022 0.066 0.145 0.276 
G 0 0.322 0.543 0.777 1.026 
H 0 0.418 0.715 1.041 1.405 
] 0.150 0.152 0.153 0.154 0.155 
K 0.110 0.1138 0.118 0.124 0.133 


The four quantifes do, a;, cr, and cg against a, with 0 
as parameter, are plotted in Figs. 3 and 4 for w = 0.35. 
Similar plots have been obtained for the other advance 
ratios uw listed in Table 1. The assumption of a con- 
stant profile drag coefficient provides approximately 
correct cg curves only in a limited range. In reality, the 
downward trend of the cg — a curves with increasing a, 
as shown in Fig. 4, is ultimately reversed. 

In Figs. 4 and 9, the cg curves are only plotted in a 
range for which a reasonable approximation is expected. 

The rotor with pitch-cone change is assumed to have 
a 6; angle of 70°. The reduction in pitch angle Aé@ pro- 
duced by an increase in coning angle Ad» is 


Aé ='‘tan 63 Ado = 2.7 Ado 
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The pitch setting is assumed to result in a blade pitch 
angle of @ = 9° at a coning angle a) = The rela- 


tion between @ and dp, therefore, reads 


= fo 
09.0 


= 


24 = 2.7ao 


Introducing this relation in the a) — a diagram, Fig. 3, 
upper portion, gives correlated values a, yu, and @ for 
the rotor with pitch-cone change and reduces the num- 
ber of independent variables from three to two. The 
three quantities a), cr, and Cg may now be determined as 
functions of a and » and may be compared to the corre- 
sponding functions for a fixed-pitch rotor, as is shown in 
Figs. 6 to 10, assuming 6 = 9° for the fixed-pitch rotor. 


LONGITUDINAL STABILITY 


Using as reference axes body axes through the center 
of gravity, the x axis pointing in the direction of undis- 
turbed flight (see Fig. 1), the equations of longitudinal 
motion of the helicopter for a disturbance from level 
flight read, if only the more important terms are re- 
tained (see reference 5, page 5), 


(W/g)% + W fadt = 0 


(W/g) (@ — aV) = 22; 


4M; + 3M: + &M, 


la = 


a is the pitching disturbance velocity, and X and 4 are 
the disturbance velocities in x and z direction. 

These three equations express the force equilibrium 
in x and gz direction and the moment equilibrium 
about the pitching axis. 

2 — &V isthe vertical acceleration in space which is 
different from Zz because of the angular velocity @ and 
the translational velocity V of the reference coordinate 
system. The three pitching moment derivatives M,, 
—M:, and —M, are velocity stability, static or attitude 
stability, and pitch damping. The derivative 
—Z:is the vertical damping. 

In a system of dimensionless quantities, using the 
rotor radius R as unit of length, the aircraft mass W/g 
as unit of mass, and the quantity mR/ V as unit of time 
(W/g)/pAR, the above equations read 


t+ (m/2)c, ff adi = 0 


Z — am = 22; 


force 


with m = 


t£M,+2M:+ aM, 


1% = 


The dimensionless quantities (symbols with a bar) 
are determined by the following equations: 


= M,/ARpV 
(M;/A RpV) 
—(M,,/AR%V) 


Velocity stability, 17, 
Static stability, —M; = 
Pitch damping, —7, = 
Vertical damping, —Z, = —(Z./ApV) 

Further, 
r = V)x 


m; x« = (m 


AERONAUTICAL 
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a= (mR/V)a; I=] (w/g)R? 

7 m2 m W m 
fe = + — sa % 
V? 2 (V?/2)Ap 2 
i = (V/mR)t 


The longitudinal motion is stable if all coefficients oj 
the frequency equation 


An‘ + Br? + OCA? + D+ ELE =0 
are positive and if 
BCD — AD? — BE>0O 


In most cases the middle term is negligible and the sta. 
bility criterion reduces to 


CD — BE>0 


The coefficients of the frequency equation read 


Az=TI; B= -M,-Z.l1 
C= M2: — Mim; D = (m/2)c,M 
E = —(m/2)c,M,Z: 


Inserting these expressions in the simplified stability 
criterion gives 


—M, — (I/m)Z2 > 0 


For zero moment of inertia of the helicopter, the longi- 
tudinal motion is stable if velocity stability, static sta 
bility, pitch damping, and vertical damping are posi- 
tive. For finite moment of inertia, the longitudinal 
motion is stable if this moment of inertia is smaller 


than the value corresponding to the equation 


I/m = 


—M./Z2 


This quantity is for the rotor with pitch-cone change 
plotted against yin Fig. 13. 


STABILITY DERIVATIVES 


During attitude or velocity changes of the helicopter 
the aerodynamic rotor torque varies. The change in 
rotor speed depends on the engine characteristics for 
Two alternative assumptions 
(B) con- 


fixed throttle position. 
will be made: (A) constant engine torque; 
stant engine r.p.m. The assumption of constant torque 
is approximately valid for reciprocating engines and for 
pressure jets at the wing tips. Actually, the torque of 
a reciprocating engine with fixed throttle position will 
slightly increase with r.p.m. because of the effect of the 
blower, and the torque of a pressure jet will slightly de- 
crease with r.p.m. In case a r.p.m. governor is con- 
nected to the power regulation, the assumption (B) is 
Actually, the rotor moment of inertia is also of 


attitude or velocity 


valid. 
some importance. For slow 
changes of the helicopter, however, the rotor moment of 


inertia may be neglected, and the constant torque 
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LIFTING ROTOR WITH 
assumption will be a good approximation. For rapid 
changes of attitude and velocity, the constant r.p.m. 
The stability derivatives are 


assumption 1s valid. 
(A) of constant 


under assumption 


first given 


torque. 

The torque coefficient q is defined by the equa- 
tion 

QO = gAR(p/2) (V + x)? 

where g. depends on » and @. Partial differentiation 

with respect to ¥ and $, together with the identity a: = 


| V,and considering x < V, 
(V 2)u; = oe IL Or 
(WV 2)us = —(1/2) (Gee! Ven) 


The pitching moment coefficient ¢y (function of a@ and 


u) is defined by the equation 
M = cyA R(p 2)( V + x)? 


and partial differentiation with respect to ¥ and 4, con- 
sidering ¢y = 0 for the steady flight state, a: = 1/1, 
%< V,and the two above equations for uw; and u:, re- 


sults in 
me . M; a 
Velocity stability, 17; = = — — Cy, 
, ; ARp|\ Fess 
Stati bilit V7 = 
static sti ty, — = = = 
itic stability ARpV 


I Vea 
at *Ma ~~ CM yu 


Ch 
Similiarly, following from 
Z = €,A(p/2)(V + x)? 


by partial differentiation and by using the above given 


expressions for a: and u:, 


hie . 5 Z: qa 
Vertical damping, —Z: = — = Cla Cl 
Ap| 2 "oa 
Finally, 
™ Ma “ee 
Pitch damping, —./7, = — ——_— “M 
AR’*p| 2 


According to Fig. 1, assuming that the thrust passes 
in the undisturbed steady flight condition through the 
c.g. of the helicopter, and approximating thrust by lift 
force, the derivatives of the pitching moment coefficient 


read 
Cc wa = C,a tall 
Cc wa = CLAi Qh 
CM = C1a;,h 
where; = h R. Substituting 
ee Se 2. — JF 2 
Q — «CQ Me, qi oo 260 be 


Qu = (2/u?) [equ — (2co/u)] 
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Cy = 2r/p*; Cra = Wra/b 
(2/u*)[cr, — (2cr/u)]) 


CLu 
and using further the approximation 
dg = —(16u/y)(R/ V) 


(see reference 5), the expressions for our four deriva- 
tives in terms of the coefficients plotted in Figs. 6 to 10 


read 
- or. 
y ° oes — Ca <CT - 
Velocity stability, M7; = — a,h 
(2¢g/u) — CQ, Be 
. ° eae = her 
Statice stability, — WM; = — — 
ye 
c ja 
Qia + Gig 
(2¢Q/ Mm) — Coy 
— : s i 
Vertical damping, —Z: = — — X 
yu? 
CQa 2¢r 
| er + - (. = 
(2¢g/m) — Coy m 
Pitch damping, —/7, = —16erh/yu 


Under the assumption (B) of constant r.p.m., one 


obtains, with uw; = uw IV and uw: = O, the following de- 


rivatives: 


- . is - Mm CQ h 
Velocity stability, 17;= — cy, =——" 
; : a 

2 Ml 

“ . eg: — l he; 
Static stability, —M; = — ~- Cy = — —> a 

: 2 pe 
Vertical damping, —-Z,= —CTo/ be 


Pitch damping, unchanged 


A comparison of the three derivatives M,, M:, and 
Z: for the rotor with pitch-cone change and for the 
fixed-pitch rotor, for constant torque, and for constant 
rotor speed is given in Figs. 11 and 12. The pitch 
damping is not plotted because it is almost the same for 
both rotor types. 

The flight conditions assumed for the curves of Figs. 
11 and 12 are defined by the assumption of a rotor 
torque equal to the torque in hovering and of a rotor 
thrust equal to the thrust in hovering. In terms of 


coefficients this means that 
Ca/Cr = Can/Crn 


where Cg, and cr, are the torque and thrust coefficients 
in hovering. 

These conditions are indicated in Figs. 6 to 10 by 
dashed lines. For appropriate values of parasite drag 
of the helicopter, these conditions are also level flight 
conditions. The main parameters for these conditions 
may be taken from Tables 2 and 3. 
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TABLE 2 


Constant Pitch Rotor 


M a . 6 Cr 103 Ce 104 ay 
0.15 —15.4 9 4.05 0.290 2.8 
0.25 —12.6 9 3.70 0. 265 4.6 
0.35 —11.1 9 3.35 0.240 6.1 
0.45 —10.3 Q 3.00 0.215 i 

TABLE 3 
Rotor with Pitch-Cone Change 

Mu a 6 cr 10% ce 10% ay 
0.15 —13.6 8.8° 4.10 0.275 2.7 
0.25 —12.0 9.1 4.00 0.268 1.6 
0.35 —11.1 9.5 3.85 0.258 6.4 
0.45 —11.0 9.9 3.50 0.235 8.4 
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The Hot-Wire Anemometer 


in Supersonic 


Flow’ 


LESLIE S. G. KOVASZNAYt 
The Johns Hopkins University 


SUMMARY 


The hot-wire anemometer has proved to be a most valuable 

tool for turbulence research in the low-speed range. As interest 
grows in the problem of turbulence at supersonic velocities, it 
would appear natural to see whether the hot-wire anemometer 
could not also be used at higher speeds. However, a number of 
fundamental problems must be solved before the hot-wire can be 
considered as a successful instrument in the supersonic labora- 
tory. " 
The first problem to be solved is that of obtaining adequate 
frequency response of the hot-wire equipment in order to resolve 
the turbulent fluctuations in space which are swept by the wire 
at high velocity in a supersonic flow. Encouraging progress 
has been made in using extremely fine wires and in developing 
electronic equipment that will compensate for the thermal lag 
of the wire at very high frequencies without interference from 
noise 

As yet, there is only a small amount of information available 
that can be used to predict the response of a hot-wire anemom- 
eter at such speeds. Because of the great difficulties of the 
theory of real compressible flows, most of the data will have to 
be obtained experimentally. The author has conducted system- 
atic experiments to determine the law governing the heat loss 
from a hot-wire in a supersonic flow up to Mach Numbers of 2.0. 
It was found that the heat Joss, expressed in the form of a Nusselt 
Number, is a unique function of the square root of the Reynolds 
Number when free-stream velocity and density are used and the 
viscosity and conductivity corresponds to stagnation condi- 
tions. This law is valid only in the case of small temperature 


loadings. For large differences, a nonlinearity 


appears in sharp contrast to King’s law, which has been found 


temperature 


to be linear with temperature difference over a wide range of tem- 
perature loadings 

The third problem that has been encountered is that of inter 
preting the hot-wire fluctuations when the number of flow param- 
eters is large. Velocity, density, and temperature all contribute 
to the hot-wire reading, and it is believed that in some cases 
these effects may be separated. When the velocity field is a low- 
intensity turbulence (not random sound waves), the separation 


of measured fluctuations is accomplished. 


Presented at the Aerodynamics Session, Eighteenth Annual 
Meeting, I.A.S., New York, January 23-26, 1950. 

* The present work represents an integral part of our research 
program on methods for measuring turbulence in compressible 
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Navy Bureau of Ordnance Contract NOrd-8036, which is under 
the technical supervision of the Applied Physics Laboratory. 


flows. 
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bauer, F. H. Clauser, and to S. Tornmarck, P. Iribe, and Miss P. 
Clarken, who helped with the experimental, computational, and 
photographic work. 
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NOTATION 


ny = wave length of a sine wave 
= frequency (cycles per sec.) 
U = velocity 
Uo = characteristic velocity of the wire 
g = temperature (absolute) 
AT = T,— 7, 
Rw = resistance of the wire at operating temperature 
k = heat conductivity 
p = density 
Ap = density fluctuation 
= viscosity 
Mors = viscosity at 273°K. 
p = pressure 
Bi = heating current 
W = heat input into the wire per unit time 
H = heat loss of the wire per unit time 
Mo = thermal lag time constant 
K = thermal capacity of wire 
Nu = Nusselt Number 
l = length of the wire 
d = diameter of the wire 
M = Mach Number 
Re = Reynolds Number 
Py = Prandt! Number 
T = temperature loading factor 
C, = specific heat at constant pressure 
a = temperature coefficient of resistivity 
4 = € To 
Y = C,,/C, ratio of specific heats 
Ae = voltage fluctuation across the wire 
Am = A(pl’) mass flow fluctuation 
n = 7./To 
1, B, ¢ / 
Ae ae fy = dimensionless constants 
A", B” \ 
F = dimensionless mass flow fluctuation sensi- 
tivity of the wire 
G = dimensionless stagnation temperature fluctua- 
tion sensitivity of the wire 
dy = {1+ a(Ty —- 273) | = (R, — Ro ;)/Ro: 
Aw’ = (Rw — R.)/R. 
R>; = wire resistance at 273°K. 
a, = 1+ al(7T. — 273) = (R. — Rey) /Ra 
Re* = B2/A? 
y = (V Re + 2V Re*)/(V Re + V Re* 
Z = 1/(1 — V Re*/Re) 
P = 1/({1 + (vy — 1)M?] 
1+ [(y — 1)/2]M 
O = 
‘ 1+ (y — 1)M 
Subscripts 
e = equilibrium temperature (unheated) 
w = wire temperature 
0 = stagnation temperature 
1 = free-stream condition 
2 = condition behind normal shock 


= mean 
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WIRE DIAMETER = .0OOI5" = 3.8microns 
OR APPROXIMATELY 7X WAVELENGTH 


OF YELLOW LIGHT 


M=I14 
R=30 


100 X MEAN FREE 
PATH 








THE HOT WIRE AND THE FLOW AROUND 
IT IS STILL VERY SMALL COMPARED TO THE 
SCALE OF THE TURBULENCE TO BE MEASURED 


Fic. 1. 


(I) INTRODUCTION 


HE HOT-WIRE ANEMOMETER has become the stand 
a tool of turbulent research in low-speed flow. 
When high-speed and supersonic wind tunnels appeared, 
the interest was focused more on the effects of com- 
pressibility than viscosity. the 
cumulation of a wealth of data on supersonic flow for 


This has led to ac- 
which no quantitative measurements were made of the 
effects of the viscosity and, in particular, of the proper- 
ties of the turbulence that was present. 

The natural development now calls for information 
on turbulence measurements in supersonic wind tun- 
nels just as was the case with low-speed wind tunnels 
in the last decade. 

The feasibility of using hot-wires in a supersonic 
flow was first demonstrated by Dryden and Schubauer 
in 1946, when they operated a 0.0003-in. diameter 
tungsten wire in the Aberdeen wind tunnel and ob- 
served fluctuations withit. (This work is unpublished. ) 

A closer examination of the problem made it appar- 
ent that there are three important problems to be 
solved before actual turbulence measurements can be 
obtained in supersonic flow: 

The first: To extend the response of the hot-wire 
and its associated equipment of the much greater fre- 
quencies of supersonic flow. 

The second: To determine the law for the loss of heat 
from the wire that will replace King’s formula used in 
incompressible flow. 

The third: 


in a compressible flow where there are three parameters 


To interpret the measurements obtained 


of the flow instead of the velocity alone, as in the low- 


speed case. 


SCIENCES—SEPTEMBER, 1950 
Before embarking on these three important prob- 
lems, it is worth while to obtain an idea of the order of 
magnitude of the difficulties involved. 
Fig. | shows an approximate picture of a hot-wire 
1.4 and 
Since 
the drag of the wire is rather large, the detached shock 


with the diameter equaling 0.00015 in. at \/ 
operating at a Reynolds Number of Re = 30. 


wave is nearly normal for that portion of the flow which 
passes near the wire. The mean free path is aboy 
one-twentieth of the diameter of the wire, and the 
shock wave thickness is of the order of one-third to 
one-fourth of the wire diameter. If we consider the 
wire plus the perturbed flow as an integral unit, even 
the fine structure turbulence will be so large compared 
to this dimension that for all intents and purposes jt 
can be considered simply as a change in the main flow. 
The characteristic time or frequency of the flow around 
the wire as estimated from the incompressible von Kar. 
man vortex street frequency, f ~ 0.2 L’ d, is of the order 
of 10 megacycles. This demonstrates that the transient 
response of the hot-wire is dominated by the standard 
heat lag and not by the time for the readjustment of the 
flow. 
(Il) FREQUENCY RESPONSE OF HoT-WIRE AND 
EQUIPMENT 


The frequency response requirements are intimately 
connected with the spatial resolution expected from 
the turbulence measuring instrument. Low intensity 
turbulence can be visualized as a slowly varying pat 
tern carried by the main flow, and the variations in 
space are converted to variations in time. A_ sine 
wave pattern with wave length \ will appear as an os- 
cillation with frequency f = l/r. This means that if 
the spatial resolution is roughly the same for subsonic 
and supersonic flows, the corresponding frequency re 
quirements increase proportionately to the mean ve- 
locity. 

If we choose 1/ = 2 and the limiting frequency in 
the equipment is finns. = 10,000 cycles per sec.,the 
corresponding wave length is approximately 2 in. By 
extending the frequency response to 70,000 cycles per 
sec., the resolution improves, since the equivalent wave 
length reduces to 0.3 in. If we define resolution as half 
the wave length, this will mean | in. for 10,000 cycles 
per sec. or 0.15 in. for 70,000 cycles per sec. at JJ = 2 
and at conventional stagnation temperatures. 

This reasoning convinced us that frequency response 
is at a premium for measurements at supersonic speeds. 
Unfortunately, the frequency response of the wires 
It starts to fall off at sev- 
eral hundred cycles and, at the higher frequencies, 


themselves is rather poor. 


falls off inversely proportional to the frequency and time 
constant. This necessitates the use of some electronic 
compensation to restore the high-frequency response. 
However, the amount of compensation that can be used 
is limited, because the signal-to-noise ratio deteriorates 
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HOT-WIRE 


proportional to the 3/2 power of the upper frequency 


jimit. The time constant of the wires improves by 
taking smaller and smaller wires, but for this, too, 
there is a practical lower limit. 

The foregoing indicated to us the necessity of de- 
veloping (1) improved techniques for using finer wires 
and (2) electronic equipment with greater frequency 
range and lower noise level. 

With regard to the use of finer wires, the stress in 
such wires is proportional to stagnation pressure and 
to the length-diameter ratio. The hot-wire method 
requires a large temperature coefficient of resistivity ; 
therefore, preference must be given to pure metals. 
An extensive inquiry was made about obtaining wires in 
small diameters. The best wire we were able to pro 
cure, handle, and operate in supersonic flow was a 
(.00015-in. tungsten wire. Its diameter is approxi- 
mately 7 wave lengths of yellow light, and we have 
taken an electron microscope picture to determine the 
surface conditions. (The picture was taken at The 
Johns Hopkins University Department of Medicine, 
courtesy of T. Murphy.) The pictures show three 
samples in Fig. 2. We also used extensively a wire 
double of this diameter (0.0003 in.). Platinum wires 
proved to be too soft. 

The wires were mounted on a wedge-shaped holder, 
as shown in Fig. 3. The tungsten wires are electro- 
plated with copper up to approximately double diam- 
eter and soft-soldered to the needle tips as described by 
Schubauer Klebanoff.* The 
portion of the wire was shorter than the span between 


and uncoated sensitive 
the needles in order to avoid the influence of the curved 
shocks from the needle tips. 

The time constant of the 0.00015-in. tungsten wires 
is approximately 0.3—0.6 XK 10~* sec., depending on 
the operating condition. 

Theoretically, it can be shown that the thermal lag is 
not dependent on the form of heat-loss law. If the tem- 
perature of the wire is 7’ and the heat generated per 
unit time in the wire is W(7’, /), with J being the heat- 
ing current, and if the heat loss of the wire depends on 
its temperature and velocity // = //(T, U), the thermal 
lag for small fluctuation is characterized by a time con- 
stant M,) that can be expressed as 


Wo = K/[(Ol1/0T) — (OW/O0T)| 


where K is the thermal capacity of the wire. As long as 
the relative temperature distribution of the wire is the 
same, the transient response of the wires for velocity 
fluctuation is the same as for heating current fluctua- 
tions. This fact justifies the calibration of thermal lag 
and compensation by the square wave technique.! 
There is a small discrepancy, however, if the transient 
readjustment of temperature distribution along the 
span of the wire is taken into account, as is shown by 
Betchov.” 

The development of hot-wire equipment capable of 
compensating the thermal lag of thin wires has been a 
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cooperative project between The Johns Hopkins Uni- 
versity and the National Bureau of Standards. Fig. 4 
shows the frequency response of this equipment when 
used under actual conditions of compensation. Fig. 5 
shows the corresponding square wave response of the 
same equipment. In both cases the time constant of 
the compensation corresponds to that encountered in 
supersonic flows with available wires. The transient 
response time is of the order of 8 microsec. correspond- 
ing to approximately 0.1-in. resolution at sonic speed. 
The noise level at this wide band operation is still not 
prohibitive and is of the order of 0.1 millivolt or only 
some 50 per cent higher than the theoretical thermal] 
noise of the amplifier. A closer investigation in terms 
of the theory of information‘ reveals that we are close 
to the point of diminishing returns, and therefore it was 
felt that we are rather close to the best compromise as 
far as low-intensity turbulence measurements are con- 


cerned. 
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gree H = (T, — T.) (k + Vv: 2rkC,pd U) (1) 
a 200 eps “ta be td TAG oe . . . : 
and in the nondimensional form can be written 
a a Nu = V (2 r)Pr V Re + (1/n) (1; 
— He om LLL te on 
or for air at normal temperature 
te ree eo Ck (ee Nu = 0.690 V Re + 0.318 (1b) 
it is used for calibration purposes as f 
— — —_—_ — - t 
is : eT ; 
SL tome VSL x00 (A). He = (Te — T) (A'VU + B’) (1 
The relation ing Eq. (16); ‘is shown in Fig. 7 as King’s 
a WIV 3000 8 formula. If, in King’s formula, we equate the heat loss 
to the heat input 
oo — ™ “ ? 40a) TT P 
ST 3000 SVS \SV\/\,_ 000 RJ? = (Tw — T,) (A'VU + B’) 
Fic. 5. Square wave response of the compensation. Time con- and differentiate logarithmically, we obtain the voltage 


stant = 0.4 millisec. 


(III) Heat Loss at SUPERSONIC SPEED 


The heat-loss problem of hot-wires at supersonic 
speed has not yet been solved by theoretical means. 
The only published estimate of the heat loss is made 
for a thin strip parallel to the flow without a shock 
wave.’ Because the wire is a blunt object with a de- 
tached shock wave in which viscous effects dominate 
and for which slip-flow phenomena are present, there is 
little hope that such simplified solutions will give re- 
liable predictions of the heat losses from wires. 

From dimensional reasoning we know that the non- 
the Nusselt Number, has the 


following functional form: 


dimensional heat loss, 


Nu = Nu 


¢ 


l 
F M, Re, Pr, r+ 


where / is length of the wire, d is diameter of the wire, 
M is Mach Number, Re is Reynolds Number, Pr is 
Prandtl Number, 7 is temperature loading factor 7 = 
(Ty — T.)/ To, T is wire temperature, and 7, is equili- 
brium temperature when the wire is unheated. 

If the Mach Number and the temperature differences 
are small, if the wire is long, and if Pr = constant, this 
reduces to the incompressible case where King’s formula 


applies. King’s formula is 
40 T T T | T T T oe ee a ia 
—— ae — OH a o-9- 0015” Dia. 
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fluctuation across the wire in terms of the velocity 


fluctuations 
Ae R, — R, l AU 
é 2R. 1+ V0,/0 0 -) 
where 
Ae = voltage fluctuation 
R, = resistance of the wire at operating temperature 
R, = resistance of the wire unheated 
Uy) = characteristic velocity, V ‘Uy = B’/A’ 


It will be noticed that the sensitivity is proportional 





to the difference between wire and air temperatures 
and to a calibration factor that diminishes with in- 
creasing velocity. Since UL really corresponds to a 
characteristic Reynolds Number (of the order of unity), 
its value can be estimated with acceptable accuracy so 
the wire is not of paramount 
the heat 


that the calibration of 
importance so long as the algebraic form of 
loss is known and we are interested only in the fluctua- 
tion sensitivity. 

Systematic measurements were carried out in a 
small specially built supersonic tunnel by varying the 
ambient pressure in the ratio of 1:4, the Mach Number 
from 1.1 to 2.0, and the temperature of wire from 
equilibrium temperature (around room temperature) 
up to 300°C. 

In the supersonic flow we still can assume that the 
Prandtl Number is constant and that the geometrical 
parameter //d is large so that the end effect can be 
computed and corrected for by the simplified theory.’ 





— 


But we have to deal with a new phenomenon. The 
“air temperature’ is no longer a unique value as it 
incompressible flow. Instead, we have a range ol | 


temperatures from the stagnation temperature to the 


stream temperatures ahead of and behind the normal 
shock. 

An unheated hot-wire placed in a supersonic aif 
stream will settle down to a temperature that we call 


“equilibrium temperature.’’ Measurements were made 
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HEAT LOSS OF WIRE IN SUPERSONIC FLOW 
Fic. 7. 


of these equilibrium temperatures, and it was found 
The 
ratio of equilibrium temperatures to stagnation tem- 
perature is almost constant, independent of Mach 
Number and of Reynolds Number until the latter be- 
comes so low that slip-flow effects appear. Fig. 6 
shows the equilibrium temperature for a 0.0003- and 
0.00015-in. wire. The large wire did not show any 
Reynolds Number effect, but the smaller wire showed 
a small increase of equilibrium temperature at the 
lowest Reynolds Number. This does not account for 
the difference in equilibrium temperature of the two 
wires, which is probably caused by the fact that they 
operate in the slip-flow region where their surface con- 
dition affects the equilibrium temperature. The thin- 
ner wire (0.00015 in.) was gold plated. The difference 
in equilibrium temperatures is not caused by a strain- 
gage effect, since it is in the wrong direction. 

After information was obtained on equilibrium tem- 


that they are close to stagnation temperature. 


peratures, experiments were made to determine the 
transfer rate when the wires were heated. 

The data for small temperatute loading (7 
appropriate wire length correction can be represented as 


Nu(Re, MM) 


> 0) with 


Nu = 


The wire length correction amounted to 7-10 per cent 
in Nu, since the length-diameter ratio was 200-400. 
The Nusselt Number and Reynolds Number can be 
defined in many different ways, depending on the par- 
ticular flow properties used. The velocity may be 
evaluated in the free stream or behind a normal shock, 
since the flow around the wire is mostly behind an al- 


most normal shock. The density, viscosity, and heat 


conductivity all can be evaluated at various conditions 
varying from stagnation conditions to free-stream con- 
ditions or conditions behind the normal shock. The 
temperatures involved have all the foregoing variants 
plus the fact that we can also use the surface tempera- 
tures of the wires. 

After trying various groupings, it was found that all 
the data would collapse into a single line independent 
of Mach Number if the flow parameters were evaluated 


behind the normal shock. 


H W Ueed 
_ A " 2p2e 4. 


2 — — B” 
nw AT ko Me 


(3) 


(p; and U, could equally well have been used because 
their product is equal to p.l/2). Since the heat con- 
duction k and viscosity u vary slowly with the absolute 
temperature, 


k/Ro = p/po = (T/T )*-*® (4) 
Thus, it makes little difference whether we evaluate k 
and yw at stagnation temperature 7) or at 7», since 
these two temperatures are nearly equal. Plotting the 
data in this form, we obtained a single line independent 
of Mach Number in the range explored. 

Fig. 7 shows the data with the straight line and its 
standard deviation. The only remarkable difference 
between low-speed and supersonic heat-loss curves is 
that the intersection with the Ni» axis changes from a 
positive intercept at low speeds to a negative intercept 
for supersonic speeds. Since negative Nusselt Num- 
bers are implausible, this indicates that the whole em- 
pirical law breaks down when it is carried to Reynolds 
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71 It is possible that the larger temperature loading jn. 
a [7 fluences the position of the shock wave and, therefor 
7 ’ 
7 the total entropy change through the shock. Hoy. 
9 ot — ever, the results, and especially their independence oj 
Pa | Mach Number, are not fully understood. 
» 4-4 “ ' 
#\ oa With the knowledge of these experimental results 
7 1 loa _| we can state the empirical form of the heat loss at super. 
ar I per 
. ©. 4 FA er aM sonic velocities. 
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2 1 If we equate the heat loss with the heat generated 
TA W = R,J? and take the logarithmic derivative of the 
Ht rat ee equation, we can determine the sensitivity of the wire 
a | |r 8 | | to fluctuations. Since p and UL’ appear only as the 
) I t &s &£ Ss s&s FF £ + Bw Ww 


HEAT LOSS vs. TEMPERATURE LOADING 
FIG. 8. 


Numbers approximately equal to unity. In this case, 
the mean free path becomes nearly equal to the diam- 
eter of the wire, and we can no longer expect the same 
type of law to hold. 

An explanation of the empirical law that has been 
obtained will not be attempted, but a “‘rationaliza- 
tion’’ can be offered as follows: The Mach Number 
behind the normal shock does not vary greatly. The 
heat loss is proportional to the square root of mass 
flow just as in incompressible flow, and, since the heat 
transfer takes place mainly on the blunt front part of 
the cylinder where almost stagnation conditions pre- 
vail, it is reasonable to evaluate heat conduction and 


viscosity under stagnation conditions. The empirical 


law 
HI V(Up)d 
Nuoy = —— =f E Bs. < (5) 
wl AT Ro Ko 
A = 0.580 and B = —0.795, with a standard deviation 


of the order of 2 per cent. 

If the temperature loading of the wire is increased, 
the Nusselt Number decreases contrary to the incom- 
pressible case, where King found linearity of the heat 
loss with temperature difference over a great tempera- 
ture range (up to 1,000°C.). The observed effect is 
a decrease of Nusselt Number with increased wire tem- 
perature contrary to the trend that would be expected 
either from radiation or from increased heat conduc- 
On the 
other hand, no systematic variation of this effect was 
detected with Mach Number or Reynolds Number. 
Fig. 8 shows the results obtained at various Mach Num- 
bers all fitting well the curve 


tivity at higher wire surface temperatures. 


H(r) = (OH/Or), «9 t(1 — Cr) (6) 


with C = 0.18. 


product m = pl, we cannot separate the effect of dens. 
ity and velocity fluctuation. As shown in Fig. 6, the 
ratio 7/7) = 7 is practically constant; consequently, 
fluctuations in equilibrium temperature and in stagna- 
tion temperature are indistinguishable. 

Omitting the algebraic steps, we obtain, for slow fluc. 
tuations (or for compensated wire in general), 


Ae/é = —F(Am/m) + G( AT /T») (8) 
where 
Ae = voltage fluctuation across the wire 
é = mean (d.c.) voltage drop across 
the wire 
Am = A(pl’) = mass flow fluctuation 
m = pU = mean mass flow 
AT) = stagnation temperature fluctuation 
T, = mean stagnation temperature 
F/G = positive constants depending on 
the operating conditions of the 
wire 
F l rnd 
1 — C’a,’(1 + ay)/g 2 (Sa 
G 
a 
1 — C’a,’(1 + ay) /g 
ng ; pF 
( — 0.384 y ay’ — nC’ax (Sb) 
l1+a, 
ze Cc ; R, — R, 
C= = ia = 
1—C R, 
where 
Ry — Reis R, — Ros 
cy. = » & > 
Ros Ros 
. = aT, 7 = ‘~ To 


Roz3 is the resistance of the wire at 273°K., where a is 
defined 


R = Rozs [1 + a(T ao 273)| 
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HOT-WIRE 


the only data entering from calibration are in z and y. 

When the heat loss is expressed for small 7 [Eq. (5)], 
Nu = AV Re + B 

If the constant Re* = —B? A? is introduced, it be- 

comes 


Nu, = A[V Re — V Re*| 


expressing the fact that the straight line in Fig. 8 does 
not intersect at zero; the coefficients y and z then be- 


come 
y oe V Re — 2V Re* 
"1 — VRe*/Re V Re — V Re* 
If the intercept in the Nu» = Nu(Re) plot vanishes, 
Re* becomes zero and z = | and y = |. Since the heat 


loss through the ends lifts the whole curve, it is possible 
to use the diameter-length ratio that gives the proper 
correction to lift the heat-loss curve the right amount. 
This is impossible for the low-speed case, since there 
the intercept is already positive. 

As we see from Eqs. (Sa) and (Sb), F increases with 
increasing temperature of the wire and G decreases 
slightly. This corresponds to similar behavior at low 
speed when, with decreasing wire temperature, the 
response of the hot-wire becomes more that of a resist- 
ance thermometer than an anemometer, as shown by 
Corrsin.’ Fig. 9 shows a practical example of a tung- 
sten wire with stagnation temperature 300°K. and & = 
0.0035 per °C. Here we have used values for y and z 
corresponding to realistic operating conditions. We 
see that, by varying the wire temperature, we are 
able to vary the relative sensitivity for temperature 


and mass flow fluctuations. The ratio # G varies al 


most an order of magnitude in the useful range. If we 
measure the mean square fluctuation 
Aé\’ | Am\  ( ATo\ Am AT, 
= FF + G* 7 — 2ru a 
é m To m Iv 
(9) 


the three unknowns ( Am)?, (AT))?2, and Am AT) can be 
determined from three measurements at different wire 
can be obtained by 


temperatures. More 


taking a larger number of readings and using a least- 


accuracy 


square method. 

If we want to measure the cross component of the 
turbulent velocity, we may use an X wire probe and 
calibrate it just as at low speed, since its calibration is 
essentially a geometrical calibration. However, care 
should be taken when the individual wires are along 
Mach waves. 


(IV) INTERPRETATION OF THE MEASUREMENTS 


As we saw above, it is possible to determine the mass 
flow fluctuation, stagnation temperature fluctuation, 
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and the correlation between the two by varying the 
temperature of a single hot-wire. 

In general, the flow conditions at a point in a com- 
pressible flow are characterized by three variables. 
One simplifying assumption will reduce it to two (e.g., 
The hot-wire gives information only 
mass flow and stagnation tem- 


isentropic flow). 
about two parameters 
perature. If we write the definition of the stagnation 
temperature 


C,T = (U? 2) (10) 


+ ly (vy — DI(p/p) 


logarithmic differentiation, we obtain an 


Introducing the Mach 


then, by 
equation for the fluctuations. 
Number and using the gas law, we obtain 


y=) eR AU 
1+ a Me) —— = (y¥ — 1) 


lo l 
A A 
eens 
pb 3B 


Eq. (11) connects the fluctuation in the four variables 
Ty, U, p, and p. If we know three of them, we can de- 
termine the fourth. Since the hot-wire can give only 
two fluctuations, A7) 7) and Am m, we need a simpli- 
fying assumption that establishes one link between the 
four variables, or we need an independent measurement 
giving a different kind of information. The isentropic 
law would give 

Ap/p = y(Ap/p) (12) 
but it is unlikely that turbulence can be approximated 
by an isentropic relation. 

In fluctuations obtained in a wind tunnel, we expect 
two distinctly different types of fluctuations. One is a 
field of random sound waves that is practically irro- 
tational, and the other is a typical incompressible but 
rotational turbulence flow field. 

For true sound waves, the isentropic assumption 


seems to be reasonable. If the intensity of turbu'ence 
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Splitting the mass-flow fluctuations into density and velocity 
fluctuctions by assuming that pressure fluctuations con be neglected. 


Fic. 10. 

is low (order of | per cent or less), the velocity fluctua- 
tions themselves have a Mach Number of the order 
of 0.01—-0.02 even for moderate supersonic mean-flow 
Mach Numbers. At such low Mach Numbers, the 
turbulence pattern is just carried along by the super- 
sonic flow and in itself is close to an incompressible 
turbulent pattern. In this case, the pressure fluctua- 
tions are of the order of Ap = (1/2)p(AU)?, and we 
can reasonably assume that the relative pressure fluc- 


tuation 
Ap/p < AU/U 


In this case, we can neglect the pressure fluctuations in 
Eq. (11), and we obtain the density and velocity fluc- 
tuations separately. The stagnation temperature may 
still fluctuate, because the turbulent pattern carries 
along a temperature “‘spottiness’’ because of its past 
history of dissipation. 

Therefore, if the velocity fluctuation is a small per- 
centage of the main flow velocity, if the Mach Number 
is moderate, and if the fluctuations are due to true 


turbulence, we should use 


AU/U = P(Am/m) + Q(AT)/To) (13) 

and 
Ap/p = (1 — P) (Am/m) — Q(AT)/To) (14) 

with 
P = 1/{1 + (vy — 1)M?] (15) 

and 
QO = {14+ [(7 — 1)/2) 473} /[1 + (vy - 1) .M*| (16) 


Fig. 10 shows P and Q as functions of Mach Number. 
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It is remarkable to note that Eq. (8) has similar for, 
to Eq. (14), which indicates that, if we chose the oper. 
ating temperature of the wire for a particular Mach 
Number so that 


F/G = (1 — P)/Q 


then 


Ap/B = —(Q/G) (Ae/é) 


This indicates that the voltage fluctuation across the 
wire is sensitive only to the density fluctuation. Be. 
cause of the difference in sign, there is no operating 
condition at which the signal is proportional only to the 
velocity fluctuations. On the other hand, it is possible 
to obtain the velocity fluctuations with the combined 
output of two wires operated at different temperatures 

If we have determined the mean square fluctuations 
both in mass flow and stagnation temperature and also 
the correlation between them, we are able to determine 
the r.m.s. velocity fluctuation, the r.m.s. density fluctua. 
tion, and also the density-velocity correlation. 


9 ° 


AU Am \" AT, \" T. 
oP) a ois) +o 
m Ty m Ty 

(17 

Ap\* ‘ AT, \’ 
hw G~ Pyle) + ei =) - 

p m L 

1 — PO ear (18 
m 1, 
ApAU Am AT, \’ 
——— = P(l—P)(—-) -@( =") + 
pl m T 
a — 2P)Q AMAT Ug 
mT, 
CONCLUSIONS 
(1) It is empirically established that the heat loss 


of the hot-wire at supersonic speed depends on the mass 
flow and stagnation temperature heat conduction and 
viscosity, and it is nonlinear with temperature loading 

2) The frequency response of wire plus electronic 
equipment can be improved to give a space resolution 
of the order 0.1—0.2 in. at moderate supersonic veloct- 
ties. 

(3) The r.m.s. 
stagnation temperature fluctuations and the correla- 
tion between the two can be determined by taking sev- 


mass flow fluctuations and r.m.s. 


eral readings with a single wire only by varying the 
operating temperature. 

(4) The data obtained with the hot-wire can be 
fully interpreted if we have some restricting assump- 
tion between the four pertinent variables 7), p, U, and 
p. If the turbulence level is low and the pressure 
fluctuations can be neglected, the evaluation gives 


(Continued on page 584) 
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Characteristics of the Turbulent Boundary 
Laver on a Flat Plate in Compressible Flow 
_from Measurements of Friction in Pipes’ 


HANS ULRICH ECKERT?t 
Air Materiel Command 


(1) ABSTRACT 


With the development of high-speed aerodynamics, the in- 
fluence of compressibility on the boundary layer becomes of 
growing importance. In the case of laminar flow, this influence 
has been treated in numerous papers, while for the turbulent 
boundary layer, which is less convenient to treat analytically, the 
literature is rather scarce. From the practical standpoint, how- 
ever, the latter is in all probability the more important, since the 
transonic and supersonic velocities achieved at the present time 
are, in general, connected with Reynolds Numbers that, accord- 
ing to experience with incompressible flow, should lie well in the 
turbulent range. In cases of extremely high supersonic speeds 
only, since they will be restricted to gases of low density, Rey- 
nolds Numbers will be low enough so that extensive laminar boun- 
dary layers can be expected. Hence, extensive turbulent boun- 
dary layers will mainly occur in the low and medium Mach 
Number ranges. Therefore, it seems to be worth while to ap- 
proach the problem by starting with formulas for incompressible 
flow and to allow for the heat due to friction by introducing 
simplified thermodynamic Attempts of this kind 
have been made by von Karman, Frankl, and Tetervin. They 
require, however, certain arbitrary assumptions regarding ap- 
propriate values for density and viscosity (values at the wall sur- 
face, in the free stream, or some mean values) which can be estab- 
lished only by experiments. The paper attempts to 
eliminate the arbitrary character by applying the results of tests 
in pipes to the plate problem by a method analogous to that 
used first for the incompressible case by Prandtl and von Kar- 


relations. 


present 


man. 
(II) NOTATION 
“u = velocity inside boundary layer parallel to surface 
(representing mean value over time) 
U = velocity of free stream in plate flow or velocity at axis 


in pipe flow (representing mean value over time) 


(subscript U refers to conditions of free stream 
and axis of pipe) 
x = distance from leading edge of plate 


normal distance of flow element from surface of plate 


¥ - 
or pipe 
6 = absolute thickness of boundary layer 
6* = displacement thickness of boundary layer 
38 = momentum thickness of boundary layer 
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d = inner diameter of pipe 
n = exponent of velocity distribution power law (subscript 
n refers to corresponding velocity profile ) 
Pp = static pressure 
p = local density 
T = absolute temperature 
R = constant of ideal gases 
uw = absolute viscosity 
w = exponent of viscosity variation power law 
Cp = specific heat at constant pressure 
oY = ¢p/cy = ratio of specific heats assumed to be 1.400 for 
air 
T = shear tension at the wall surface 
9 z 
Cc = — rdx = mean coefficient of friction at the 
pul 2x 0 
plate 
Cy’ = 2r/pyU? = local coefficient of friction 
M, = U/WV yRTvy = Mach Number of free stream 
Re = Upyx/uy = Reynolds Number of plate flow based on 
free-stream values 
A = cross section area of boundary layer on plate and in 
pipe 
eA 
i = udA = mean value of u over A 
ew & 
p = pdA =mean value of p over A (mean static density) 
74 
Pp = pudA = mean density of mass flow (mean 
aA 
eo 
flow density) 
1 7A 
p = pu’*dA = mean density of momentum flow 
ped 
(mean dynamic density 
T = p/pR = definition of mean temperature based on the 
assumption of a mean state across pipe section 
Me = value of u based on 7 
Rea = tipd/u = Reynolds Number of pipe flow based on 
mean values defined above 
r = &r/pi? = coefficient of local friction in pipes 


Subscript 7 refers to incompressible flow 
Subscript 0 refers to stagnation conditions 


4 = pP/Pt 

y = y/é 

V = 4u/U 

A = 6*/6 

6 = 0/65 

K = [(7 — 1)/2])M,? 

a = (T, — Tv) To=2 K (+k 

K = 0.4 = constant of mixing length theory 

s = 11.5 | Constants of laminar sublayer for incompressible 
Ff = (0.289 ( flow 


ams kV 2 Cyx’ = KV p,l? T 
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DEFINITION AND DETERMINATION OF MEAN 


DENSITIES FOR PLATE AND PIPE FLOW 


(IIT) 


(A) Basic Assumptions 

Before discussing the main problem it seems wise to 
define the different mean values for the boundary- 
layer parameters which will be used and to derive the 
formulas to compute them for plate and pipe flow. We 
start with the following assumptions: 

(1) No heat is exchanged between the fluid and the 
wall, and exchanges of heat and kinetic energy between 
the elements of the fluid balance each other. Then 
the total energy per unit mass is constant and equal 
to the value in the free stream and in the axis of the 


pipe. 
(u?/2) + cD = (U?/2) + GT, (1) 


(2) The velocity distribution in the boundary layer 
is given by a power law, 


u/U = (y/6)'" (2) 


with » constant. For fully developed pipe flow, 6 is 


replaced by d/2. Deviations from the given profile 
near the wall (laminar sublayer) and near the free 
stream are neglected. 

(3) The static pressure across the boundary layer 
is constant and equal to its value in the free stream 
and in the pipe axis. 

p = pu (3) 
For plate flow py is constant. 

From Eq. (3) the relation between density and tem- 
perature follows for the plate flow and a given cross 
section of the pipe flow. 


p/pu = Ty/T (4) 


or, with the use of Eqs. (1) and (2), 


Introducing the Mach Number with 
U? =o,(7 — 1) Ty M,? (6) 


we obtain for Eq. (5) 


=a) > de 
a te ue|1- (2) [} (7) 
pu \ 2 6 
or, if setting for brevity, 
fs) v — 1 
cae tay: + M,? = K 
Pr 6 2 
(Va) 
l 
P aim 
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Fic. 1. Dependence of the ratio mean static density to fre 
stream density on exponents of velocity distribution law for yar 
ous Mach Numbers. 


(B) Mean Static Density 


By integrating Eq. (7a) from Y = 0 to Y = 1, we ob- 
tain for the plate flow a term that represents the mean 
density over a cross section in relation to the free 








stream density. 


plate ft dY | 
r = : = ) 
. o 1+ A(1 — Y°") 


(The subscript denotes the shape of the velocity | 


profile.) We call P the nondimensional mean static 
density and will see that the other characteristic mag- 
nitudes for plate and pipe flow can be deduced from this 
term. 

The integral of Eq. (8) has been solved analytically 
and evaluated for a wide range of K and n values. For | 
n = '/o, 1, and 2, direct solution was found feasible, 
while for even values of 7 between 2 and 22 the inverse 


a 
J YdP with Py = 1/(1 
Po je 
is equal to | -f PdY. 
0 


values of ” are found by interpolation. 


+ K) has been solved which 


Values of P for intermediate 


The results are given in graphic form in Fig. 1, with 


the Mach Number as parameter. The exponent ” = | 
is sometimes used as a rough approximation for the 
laminar boundary-layer profile at the plate, while for 
turbulent flow m may vary between 5 and 10, with” = 
7 asa most commonly used value. 

To obtain the mean static density for pipe flow we 
have to solve the integral 


pipe | ( | V)dVY 
P, = 2 : = d 
I i+ Ka — ry”) 


which can be split into 


| “I dY 
" J 1+ K(1 


Vdy 


ad | 
ee 1+ K(1 — 
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The first integral has already been solved for the plate. 


second one can be reduced to the first by placing 


The 
Carrying out the substitution and setting for 


y=Z. 


brevity here and in the following, 


1+¢K0-Y°") = F(Y) = 14+ K(l — 2’) = GZ) 


we obtain, from Eq. (9), 


pipe "1Z"-'dZ ~— (2n) {3 | 
P, = In a — - — (Ya) 
7 y G(Z) y 4 J0 G(Z) 


Evidently, the second integral denotes half the plate 
value for the case 2”. 


Thus, we have 


pipe plate plate 


| = 2 Ps 1 Po, (10) 


(C) Mean Density of Mass Flow 


Replacing the velocity distribution by the mean 
velocity and requiring continuity of mass, we obtain the 
expression for the mean density of the mass flow. 


With a/U =V for the plate, 
: 1 
PV = pp = frvay (11) 
0 
As 
Vy. = n+ 1) (12) 
we have 


n+ | 1y’/*qy 
P, = (1:3) 
0 an 


Again substituting Y = Z”, we obtain 


Z" 
P, = (2 + I (13a) 
G(Z 
Setting, temporarily, m = n’ — 1, we have 
of 1 
_—. = & — dz = P,, (13b) 
9 G(Z) 
or 
plate plate 
P, = Pi: (14) 
For the pipe, 
4 "1 y'/n(] Y)d) n 
rvs = 2 (15) 
0 F(Y) 
with 
a 2n° 
V, = (16) 
(mn + 1) (2n + 1) 
we obtain in the same way 
pipe 2n+1 plate n + 1 plate 
P, = P — P. (17) 


+1) 2n + 
n ‘ n — 


LAYE 


ur 
cr 
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(D) Mean Dynamic Density 


From the assumption of a mean momentum flow, we 
can define a mean density of momentum flow or mean 


dynamic density. For the plate, 


1 
PV? = P/V? = f Pid) (18) 
0 
with 
V,2 = n+2 (19) 
we obtain 
plate shine - 
2U) 
a ™= P. 4-2 
For the pipe, we obtain with 
2n? 
V2 = (21) 
(n + 2)(2n + 2) 
pipe 
Qn + 2 ¢ Plate n+2 ed (29) 
i ui 2 
P, = P.a2 — 2n 4-2 
n n 


(E) Displacement Thickness 


The displacement thickness for the flat plate is de- 
fined by the expression 


75 
pu 
§* = J (: = Jas (23) 
0 pul 
or in nondimensional form, with 6*/6 = A, 
ag | 
A = / (1 — PV)dY (23a) 


Jo 


12), and ( 


From Eqs. (11), ( 14) of Section (C), it follows 


that 





Fic. 2. Dependence of the ratio displacement thickness to 
full thickness on Mach Number for various exponents of the 
velocity distribution law. 
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Plate] ~ 
n 


n+ 1 Pats | 24) 


A,=1- 
Numerical values for A are given in graphical form 
vs. M,, with n as parameter, in Fig. 2. 


(F) Momentum Thickness 


The momentum thickness (more correctly momen- 
tum defect thickness) for the flat plate is defined by the 


expression 
75 
pu u 
e= / ( _ Jay (25) 
0 pul l 
or in nondimensional form, with 3/6 = 0, 
“ 
0 = / PV(1 — V)dY (25a) 
0 


Using Eqs. (11), (12), and (14) from Section (C) and 
Eqs. (18), (19), and (20) from Section (D), we obtain 


plate plate | 


n 


n 
n + 2 Pate 


“ies n 4 1 | ae ‘ (26) 


0 
Numerical values for 0 are given in graphical form vs. 
M,, with as parameter, in Fig. 3. 

It should be noted that in the same manner a dissipa- 
tion thickness can be defined and computed. 
DISCUSSION OF MEASUREMENTS ON 
COMPRESSIBLE FLOW IN PIPES 


(IV) 


The results dealt with in Part (III) merely take into 
account the density variation within the boundary layer 
and say nothing about its actual extension. They 
thus represent what may be called an interior compres- 
sibility effect. To obtain information about the influ- 
ence of Mach Number on 6 itself or an exterior com- 
pressibility effect, we either have to make additional 
and rather arbitrary assumptions or look for some ex- 
perimental data on 6* or the friction coefficient. As 
measurements on plates are not available, we will use 
data from experiments on compressible flow in pipes. 
To make them applicable to the plate problem, we have 
to develop a suitable analogy between pipe and plate 
flow which includes compressibility effects. The experi- 
ments to which we refer have been described in a paper 
by Fréssel‘ and a more recent one by Keenan and 
Neumann.’ [In both works, flow of 
atmospheric air through smooth cylindrical pipes at 


experimental 


subsonic and supersonic speeds has been investigated. 
As, in the majority of tests, the stagnation tempera- 
ture of the air did not differ much from the room tem- 
perature, the requirement that no heat exchange takes 
place between the wall and the fluid should be fairly 
well fulfilled. 

Another necessary condition, the existence of which 
cannot be confirmed so readily, is that the friction coeffi- 
cient has to be determined for fully developed turbulent 
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Fic. 3. Dependence of the ratio momentum thickness to full 


thickness on Mach Number for various exponents of the velocity 
distribution law. 


+ 





pipe flow (that is, extension of the boundary layer to the 
center of the pipe; in the following for briefness thisis 
called “‘pipe flow’). The establishment of such a flow 
requires a minimum length of the pipe which depends 
primarily on the entrance conditions. On the other 
hand, wall friction for compressible flow may lead to 
choking—that is, the flow tends to attain the state of 
Mach Number | regardless of whether the entrance 
Mach Number is subsonic or supersonic. To attain 
pipe flow at a desired local Mach Number, the Mach 
Number at the entrance has to be lower in the subsonic 
case and higher in the supersonic case than the local 
Mach Number. For entrance Mach Numbers in the 
neighborhood of 
These conditions have been investigated 


1, no pipe flow is obtained before 
choking. 
analytically for various entrance Mach Numbers and 
Reynolds Numbers by Young and Winterbottom." 
The results of these investigations indicate that at a 
Reynolds Number R,, = 2 X 10°, which is the average 
occurring in the experiments of references 4 and 5 and 
for the supersonic case that is the critical, pipe flow can 
be established for entrance Mach Numbers above 1. 


before choking occurs. For practical application, this 


result has to be used with caution, since no allowance f 


is made for a laminar part or special disturbances at 
entry. 
the general behavior of the flow. 

Inspection of Fig. 7 in Fréssels paper, which repre: 
sents the pressure variation with pipe length at super- 
sonic speeds, gives the impression that no transition 
laminar turbulent occurs and that turbulence must 


have started ta the pipe entrance. This is supported 
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TURBULENT BOUNDARY 
by preliminary tests of Keenan and Neumann with 
Laval nozzles of similar shape to those used by Fréssel. 
The authors suggest that oblique shocks originating at 
the junction between nozzle and pipe may have favored 
the development of turbulent pipe flow. 

This statement is important because Frdéssel’s data 
would otherwise be rather useless for our purpose, since 
he gives only mean values for the friction coefficient \ 
and R,, which are taken over the pipe length attained in 
each particular case. 

For their main tests, Keenan and Neumann used 
more slender nozzles that were especially shaped in 
order to avoid any shock at the pipe junction. Prob- 
ably this accounts for the difference between their re- 
sults and Fréssel’s. The local friction coefficient vs. 
pipe length, as presented in Fig. 14 of reference 5, 
shows a complicated behavior near the entrance (first 
2) diameters) which may be explained by transition 
from laminar to turbulent flow in that region. 

Therefore, it seemed logical for our purpose to use the 
values at the end of the supersonic flow (25-50 diam- 
eters), since those probably most closely represent fric- 
tion coefficients for pipe flow. They have been taken 
from a curve faired by the authors through the scattered 
points of each test series. 

Fig. 4 gives a survey of the A values vs. R,, for super- 
sonic flow. Fréssel’s values are taken from Fig. 10 of 
his paper and are brought into explicit form. The 
values of Keenan and Neumann are multiplied by 
four to compensate for the different definition of X. 
The end Mach Numbers vary between 1.2 and 1.6 only. 
On this basis, extrapolation to higher Mach Numbers is 


uncertain. The dotted curve represents the relation of 



































Fic 4. Friction coefficient in pipes at supersonic speeds. 
Experimental data: XX by Fréssel‘ are mean values over 
x/d, with x/d varying between 1.6 and 27. Mach Number at 


Pipe inlet varies between 1.7 and 3.6. © by Keenan and Neu- 


mann? are local values taken at (x/d)maz 


Ma at 
(5) (5) 
No Marntes d} maz d/ maz. 
I 2.84 38 1.20 
I] 3.87 50 1.21 
III 2.91 36.5 1.38 
IV 3.14 38 1.58 
V 2.06 25 1.3: 
Relations for incompressible flow: von Karman 


Nikuradse, Eq. (27): Blasius, Eq. (28). 
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Karman and Nikuradse established for incom- 


pressible flow.® 


von 


1/Vd = —0.8 + 2 log (RaVd) (27) 
The straight line is the Blasius law 
\ = 0.316 R,4~ (28) 


For subsonic speeds, both investigations have led to the 
result that the relation between \ and R,, can be rep- 
resented within a few per cent by Eq. (27)—4.e., no in- 
fluence of Mach Number is apparent. 

For supersonic speeds, Fréssel concluded from his 
results [which, though more scattered than the sub- 
sonic data, show no clearly deviating trend from Eq. 
(27)| that the same relation is valid, and, hence, Mach 
Number has no effect on skin friction in pipes. Keenan 
and Neumann more cautiously suggest that \ might 
reach the incompressible value at great pipe length. 
This would that the Mach Number effect is 
negligible at low supersonic Mach Numbers only. 
Indeed, the results of these authors indicate that at 
higher Mach Numbers the friction coefficient might be 
Before this in- 


mean 


lower than for incompressible flow. 
fluence can be clearly established, however, it seems 
“The 


friction coefficient for compressible pipe flow is, at 


more justified to use the conclusion of Fréssel: 


subsonic and supersonic speeds, equal to the value for 
incompressible flow of the same Reynolds Number.” 

Instead of Eq. (27), we shall use as a basis for the 
further calculations the Blasius law, Eq. (28). This 
decision seems to be tenable because of the uncertainty 
of the experimental results. It will greatly simplify 
the calculations, since it involves the assumption of a 
constant velocity profile (n = 7). Up to now, there is 
little known about the influence of Reynolds Number 
on the velocity profile of compressible flow.t 

The nondimensional terms for \ and R,, are calculated 
in both papers under the simplifying assumption that a 
uniform velocity across the pipe section exists (in refer- 
ence 5, called determination of the “‘apparent friction 
coefficient’). Consequently, the differences between 
mw and +~/u2 and between the different mean values for 
p and 7 which originate from the existence of a velocity 
profile are disregarded. While the former inaccuracy 
is of the order of 1 per cent, the latter will lead to ap- 
preciable errors at higher Mach Numbers. 

+ Fréssel’s paper includes experimental determinations of the 
velocity distribution at found that 
beyond a pipe length of 36 d the shape of the velocity profile does 


subsonic speeds. Fréssel 


not alter appreciably. This agrees with Nikuradse’s tests on 
water.® Also, the shapes of the fully developed velocity profiles 
nearly coincide with those of Nikuradse and justify the approxi 
mate presentation by a power law. Within the investigated 
range of R,, from 142,000 to 275,000, 2 varies between about 7.4 
and 9.8. According to Nikuradse’s tests, the corresponding 
range in R,, for water is 105,000—2,350,000. 
increase of » with R,, seems to be much stronger than for incom- 


Consequently, the 
pressible flow. It cannot, however, be read out of the test results 
with the necessary accuracy to make the application of a com- 
plicated method appear worth while. 
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Besides the stagnation values pp) and 7), the quanti- 
ties to be measured are the static pressure along the pipe 
p and the mass rate of flow pu, which in our notation 
is identical with p#. p and a are then determined by 
the laws of conservation of mass and energy. The 
temperature of one cross section 7 is calculated under 
the assumption of a mean state as 


T= p/ pR 


From 7, the viscosity value y is obtained by Suther- 
land's or a similar relation. * 


(29) 


APPLICATION OF THE TEST RESULTS TO 
CONDITIONS AT THE FLAT PLATE 


(V) 


(A) The Principle of Prandtl and von Karman 


In order to apply the results from references 4 and 5 
to the plate problem, we remember that the fully 
developed pipe flow means a boundary layer that 
has grown from the wall to the center and has reached 
its maximum thickness, the pipe radius. Furthermore, 
we recall the hypothesis established by Prandtl' and 
von Karman? that the shear stress exerted by the fluid 
on the wall depends on viscosity, density, and velocity 
distribution in the immediate vicinity of the wall only. 
As the thickness of this effective layer is small in com- 
parison with the pipe radius, the extrapolation from 
finite to infinite radius, which is the plate, should not 
introduce an appreciable error if the other flow condi- 
tions are equal. These require, first, equal values for 
u, p, and LU’ and identical velocity profiles (for instance, 
'/; power law) for the two cases. Then, at a position 
along the plate where the boundary-layer thickness 6 
equals d 2, the velocity distributions “(y) are identical. 
Hence, a relation between shear stress and Reynolds 
Number obtained from pipe tests can be applied to the 
plate if d/2 as characteristic length is replaced by 6. 

That this principle leads to good agreement with 
experiments on incompressible flow has been shown for 
low Reynolds Numbers and constant velocity profile 
(n = 7) by Prandtl' and von Karman’, and for high 
Reynolds Numbers and variable velocity profile by 
Schiller and Hermann.* 

We assume that this principle can be extended to 
compressible flow. In this case, u and p also depend on 
y. However, according to the assumptions of Part 
(III), the functions for the pipe and the plate at the 
considered position are similar. The 
equality at the wall will therefore be fulfilled if we have 
equality at an arbitrary wall distance. For convenience 
we will use the values at the edge of the boundary 


conditions of 


layer wy and py. 


(B) Conversion of Test Results to Conditions at Pipe 
Center 
The necessary 
Reynolds Number is given, in our case, by the Blasius 
law [Eq. (28)], with \ and R,, based on j, p, and @. 


relation between shear stress and 
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For application to the plate problem, we must replace 
these mean values by the values at the pipe center 
Bu, pu, and U, 

For n = 7, with Eq. (16), is 


u/U = 98/120 = 0.817 (16a 
independent of the Mach Number. 

The ratio p/py depends on 7/Ty only, while j/y, 
depends on 7/7, as well as on the absolute temperature 
itself. However, as the considered temperature range 
is not too far, we can use with sufficient accuracy the 


relation 
h/ee = (T/Ty)* (30 


0.8. 
From Eqs. (3) and 


and take an average value of w = Then we must 
determine the ratio 7/Ty only. 


(29), it follows that 
(31) 


Numerical values for 7/7, have been computed with 
Eqs. (17) and (31) and are presented in graphical 
form in Fig. 5 as a function of the Mach Number at the 
pipe center with m as parameter. For the friction co- 
efficient, we use the definition corresponding to Fréssel’'s 
equation, 


\ = 8r/pit? (32 


Introducing Eq. (32) into Eq. (28) and writing pad/p 
instead of R,,, we obtain, for 7, 


t = 0.0395 pit?(a/tpd) (33 


or, using center values with the aid of Eqs. (16a), 
(30), and (31), 


Ku 


1/4 T [(w+1)/4]-1 
r = 0.0277 pot ( ) (*) (34 
Uprd Ty 


(C) Full Thickness of the Boundary Layer 


We consider now the case of a flat plate parallel to 
the flow whose free-stream conditions equal the center 
conditions at a pipe section. In this case, according 
to the aforesaid, Eq. (34) represents the shear stress 
for that distance from the leading edge where 6 equals 
d/2. On the other hand, the shear stress is given by 
von Karman’s momentum theorem. As UL’ is considered 


constant for the plate, we can set 
af Hv : 
0.0233 py U? ( a ) x 
l pvd 


d ff ’ 
/ pu(U — u)dy = 
dx Jo 
— (3—w)/4 
( 1 ) (35 
Ty 


After dividing Eq. (35) by pyl®, the integral becomes 
identical with the definition of the momentum thick- 
ness #. The ratio 3/6 = O has already been deter- 
mined from Eq. (26) in Part (II). From the suppost- 
tion that the shape of the velocity profile does not alter 
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ture of center in pipes vs. Mach Number at pipe center for various 
exponents of the velocity distribution law. 
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Effect of compressibility on boundary-layer thickness 
for 1/7 power law. 


Fic. 6. 
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Fic. 7. Effect of compressibility on local and mean friction 


coefficients vs. Mach Number for !/; power law of velocity dis- 


tribution 


with the plate length, it follows that 0 is not a function 


—(3—w)/4 
) (36) 


= x, 


of x, and we can write 


16 1/4 7 
i - 0.0233 ( = ) ( 
dx Upvd, Ty 


or, after integration and with Upyx/uy 


5 F r —(3—w)/5 i 
= 0.059 87 “4 — R, (37) 
x Ty 
For incompressible flow, Eq. (37) reduces with 


T/T, = 1 and 0, = 0.0972 to the well-known expression 


0.381 R,7’ (37a) 


6,/x = 


AY 


can be given. 
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(37) by (37a), we obtain, for the 


compressibility effect on 6, 


6 (2) Fy" 
5, \O Pa 


The ratio 6/6; is plotted vs. Mach Number in Fig. 6 
taken from Figs. 3 and 


Dividing Eq. Eq. 


(38) 


with values for 0/0, and 7/7; 
5, respectively. As can be seen, 
and, for 1/, = 10, provided the underlying assump- 
the 2 would be at- 
tained. The increase is considerably that 
predicted for the laminar boundary layer by von Kar- 


there is some increase, 


tions retain their validity, value 


less than 


man and Tsien.® 


(D) Friction Coefficients 
The J/ocal friction coefficient is obtained from the 
momentum theorem as 
C,’ = 20(d6/dx) (39) 
As, from Eq. (36), it follows that dé/dx = (4,5) (6/x), 
we obtain, with Eq. (37), 
C,’ = 0.09440'"(7/Ty)~ °@-”?R,~ (40) 
This reduces, for incompressible flow, to 
Cy’ = 0.0592 R,- (40a) 
and the compressibility effect is 
4 hae Be ied 
— as me (41) 
Cri 0; Ty 
The mean friction coefficient is 
C, = 20(6/x) (42) 
or, with Eq. (37), 
C, = 0.1180'"(T/Ty)~ 9-9 R,” (43) 
which reduces to 
Cy, = 0.074 R, (43a) 
and we obtain, as for the local coefficient, 
"oe ray V/s rh (3 —w)/5 
ais: ae (44) 
C fi 0; 7; 


The ratio of the friction coefficients is plotted in 
Fig. 7. It shows a considerable decrease with Mach 
Number. For /, = 10, 
compressible value is obtained, while for laminar flow 
a decrease to about 70 per cent is predicted in reference 
S for that Mach Number. 


only 35 per cent of the in- 


(VI) PROBABLE LIMITS OF APPLICABILITY 


The results obtained in Part (V) have to be supple- 
mented by some considerations concerning their field 
of applicability. So long as they cannot be based on 
rough estimations 


it 


experimental facts, however, only 
Regarding the Mach Number range, 
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is believed by the author that, beyond Mach Numbers 
of about 3 to 4, the basic simplifications will no longer 
be tolerable and that, when reliable test results for pipe 
flow at higher Mach Numbers are available, a more 
refined calculation method should be applied. There- 
fore, the present results for higher Mach Numbers can 
be regarded as qualitative only. 

Regarding Reynolds Number, the range of applicabil- 
ity as a first approximation can be set equal to that for 
which the R,~'”° 
This law assumes existence of turbulent flow from x = 0. 
In general, however, allowance must be made for the 
laminar part whose practical upper extent 
sponds to a transition Reynolds Number of about 
3X 10°. 

Some conclusions concerning the possible variation of 
this limit with Mach Number can be drawn from the 
theoretical investigations of Lin and Lees on the sta- 
bility of the laminar boundary layer in compressible 
flow.'* 14 The investigators found that, in the case of 
the heat-insulated wall (74, = 7), the critical Rey- 
nolds Number—that is, the value of Reynolds Number 
at which potential disturbances no longer would disap- 


law for incompressible flow applies. 


corre- 


pear by damping—considerably decreases with increas- 
ing Mach Number. 
tion Reynolds Number from these data, it has been 
assumed that the ratio R R = 6, which is 


ecrit 
frequently used at incompressible flow, can be main- 


To obtain values for the transi- 


etrans 


tained. 

It should be mentioned, however, that the values for 
this lower limit obtained for the special case of zero heat 
transfer seem to be of little general importance, since 
another result of Lees’ paper’! indicates extreme de- 
pendence of R,,,,, from the wall temperature. By 
heat transfer to the wall, the laminar boundary layer 
becomes stabilized, and, at supersonic speeds, laminar 
flow can be retained even at all Reynolds Numbers if 
the wall temperature is kept only 12-15 per cent below 
Zo. 

The upper limit for the validity of the R,~’* law for 
incompressible flow is found at about R, = 107. Some 
indication about the variation of this figure with com- 
pressibility can be obtained from Eq. (37), which 
gives—after extension of the left side with U’p;, uy, the 
use of Eqs. (16), (29), and (30) and 6 replaced by d,2 
a relation between the Reynolds Numbers of the pipe 
and the plate and which can be solved for R, as 


R, = 18.60(T/Ty)?**R.q* (45) 
This, in the incompressible case, reduces to 
Re @ 1S, (45a) 


As is known, Eq. (45a) connects with fair agreement 
the ranges of validity of the R,,~‘ law for the pipe 
and the R,~ ‘law for the plate. 
that the upper limit for the Blasius law is independent 
from Mach Number, we obtain, by dividing Eq. (45) 


If we assume, first, 


by Eq. (45a), 
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Fic. 8. Estimated limits of applicability for R,~'/*-law ys 
Math Number. R.,,,, = Eq. (46); Remin. = Retrans. = 6Reeit 
Data for R,.,;, from Lees.'* 
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From Eq. (46), it follows that the upper limit will 
be extended to higher Reynolds Numbers with growing 
Taking R,,,,,. = 107, we obtain, for 
Ma = 5, Renar. © 4X 107. If the upper limit for the 
Blasius law in compressible flow is higher than at ./, = 
0, which might be indicated by the test results of 
Keenan and Neumann (see Fig. 4), the upper limit for 
the R~'* law would be further increased. 

Values obtained for the lower limit at zero heat trans- 
fer from Lees’ data and the upper limit from Eq. (46) are 
plotted in Fig. 8. It is emphasized that they are of 
qualitative character only. 


Mach Number. 


COMPARISON WITH THE RESULTS OF OTHER 


INVESTIGATORS 


(VII) 


(A) von Karman 

The first known proposal to allow for the compresst- 
bility effect in the friction of the turbulent boundary 
layer originates with von Karman.’ Proceeding from 
the general relation for incompressible flow, 


VC, = 0.242/log (R,C,) (47 


von Karman relates C; and R, to values of density and 
viscosity at the wall surface and obtains, with Tyan = 


To, 
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(0.242 E (R.C;) — w log (: 4. 


y 


ut) (48) 


Eq. (48) results in a considerably stronger decrease of 
C, with Mach Number than Eqs. (43) and (44). The 
ratio C,/Cy, from Eqs. (47) and (48) is plotted for com- 
parison in Fig. 9 for R, 10° and R, = 108. The 
curves indicate a growing influence of compressibility 


with Reynolds Number. 


(B) Frankl and Voishel 


A more detailed treatment has been given to the 
problem in two papers by Frankl and Voishel,® ' who 
try to extend the mixing-length theory to compressi 


C/ = C¢ 


Sx 


The curves C,’/C;,’ for R, 10° and R, 


as is evident from Fig. 9. 


The ratio 6/6; can also be obtained from Eq. (36) of reference 9. 


stream values for R, and C,’, reads 


| «wv(2 cy’ (1 a) s j 
E + (« l 
kf 


6 V(2/C/) (1 — a) 


r R.A 


a)ite 


The 6 values obtained from Eq. (52) are sensitive 
toward variations of the C,’ values. As these have 
been found for the Reynolds Numbers under discussion 
by interpolation and are already afflicted with some 
uncertainty, it is natural that the obtained 6 values are 
somewhat scattering. Nevertheless, it is obvious from 
the plotting of the calculated points in Fig. 10 that the 
increase with Mach Number does not go on monoto 
nously but in wavy fluctuations. This result is believed 
to be due to the complicacy of the mathematical pro 
cedure, which requires expansion of functions into 
power series several times. The average trend of the 
curves is not very different from that obtained by Eq. 
(38) of this paper. 


(C) Tetervin 


Starting with the same assumptions as have been 
made in this paper, Tetervin'' has brought data in more 
implicit form, which, besides plate flow, also consider 
problems of rotational symmetry and include the exist- 
cate of a pressure gradient. The results of pipe tests 
irom references 4 and 5 have also been considered. 
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bility and include also the case of heat transfer to the 
fluid. Because of the analytical difficulties, however, 
various physical and mathematical approximations 
have to be admitted. The velocity profile is computed 
with the assumption that the shear stress through the 
boundary layer is constant and equal to the value at 
the wall surface. The nondimensional combinations 
for velocity, wall distance, etc., are based on wall values 
for w and p. The Prandtl number in the laminar sub- 
layer and its equivalent in the turbulent layer are set 
equal to unity. For the case of zero heat transfer, 
numerical values of the local friction coefficient are given 
in reference 9 at two subsonic speeds only. To extend 
the comparison of results to higher Mach Numbers, 
values have been calculated up to J, = 4.4 from Eqs. 
(11), (12), and (13) of reference 9, which, if changed to 


the notation of this paper, are as follows: 


Z 


| a Sata 85.336 _ "e¢ 
+ =e 15.556 4+ — 38.5a dz (49) 
» Zz z 
— q)'t’ (50) 
(1 — a) (51) 


10° nearly coincide with those obtained from von Karman’s Eq. (48), 
The reason might lie in the relation to wall values common to both methods. 


This equation, if solved to 6/x and with free- 


C,'xs* ) 


They are interpreted in a way that no influence of Mach 
Number on the friction coefficient exists at all, dis 
regarding the fact that these results were obtained 
with Reynolds Numbers based on mean values while 
the formulas used by Tetervin are based on free- 
stream values. The ratio of mean Reynolds Number 
to free-stream Reynolds Number, however, depends 
on Mach Number, as has been shown in Part (V); con- 
sequently, the friction coefficient, if presented as func- 
tion of free-stream values, is also a function of Mach 
Number. The conclusion drawn by Theodorsen and 
Regier,'? from tests on dises revolving in freon gas 

that the friction coefficient does not depend on Mach 
Number!!—is in disagreement with references 4 and 
5 and, with to 
supersonic wind tunnels, no longer seems to be ten 


regard more recent experience in 


able. 


If the assumption is made that C,/C;,, = 1 for the 


ratio of boundary-layer thickness, it follows that 


0,/0 (53 


which, as shown for m = 5 and 7 in Fig. 10, results in a 








582 JOURNAL OF THE 


. 


il 


—— MA 
O77 3 4 5§ 


Fic. 9. Compressibility effect on skin friction coefficients. 

Comparison of results from different theoretical investigations: 
I, von Karman; II, Frankl and Voishel;' III, Tetervin;'! IV, 
Iq. (41) and Eq. (44) of this paper 0.8. 





qo = 


much higher increase with Mach Number than Eq. 
(38). 
(D) 


Wilson, Young, and Thompson 


After completion of the work for this paper, the 
author received notice of an experimental determina- 
tion of boundary-layer profiles on a flat plate in super- 
sonic flow 


by means of pitot tubes performed by 


Wilson, Young, and Thompson." As the performance 
conditions closely agree with the requirements of this 
paper, the tests should provide a suitable basis for 
checking the theoretical results above. 

One important outcome is that the '/; power law for 
the velocity distribution, to which the results of this 
paper after Part (IV) have been specialized, is confirmed 
with fair accuracy for plates up to Mach Number 2.2 
and Reynolds Number 1.9 X 107. 

The other quantity to be checked is the boundary- 
layer thickness. As the experimental determination of 
the absolute thickness of the boundary layer depends 
on the sensitivity of the measuring instruments, it is 
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considered more reliable to compare the Calculated 
and experimental values of the displacement thickness 
The values of 6* dealt with in reference 15 are plotted 
in Fig. 11. They can be roughly divided into ty, 
groups for average Mach Numbers 1.7 and 2.0, 7 
values for M, = 


‘ 
Ne 
2.0, which include four test SeTies 
had to be shifted by amounts x, versus smaller y { 


allow for a laminar part. The x, values have been de. 
termined for each series from the differences in Rey 
nolds Numbers between Figs. 13 and 14 of refereng 
15, since no direct data are given. 

The two continuously drawn curves give the theoreti 
cal 6* values for M/, 1.7 and 2.0. They are obtaine 
from the relation 


oe eu * 


6* = (6*/6) (6/6,) 56; 4 


with the ratios 6*/6 and 6/6; taken from Figs. 2 and{ 
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Fic. 10. Compressibility effect on boundary-layer thickness 

Comparison of results from different theoretical investigations 
I, Frankl and Voishel;" II, Tetervin;'! III, Eq. (38) of this paper 
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11. Displacement thickness vs. plate length. Compar 
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Fic. 
son with experimental data from reference 15. 


been shifted by x, to allow for laminar part: 
A M, ® x, = O 
] Mee 1.9 x, = 3.3 in. 
M, & 2.0 x, = 3.6 in 
xX MaA 2.12 x_ = 6.0 in. 
+ Maz 2.18 cz, = 6.1 in. 
Present Theory: 
I M, = 2.0 R./x = 655 X 10?/in 
I] M,= 1.7 R./x = 575 X 103/in 
II] M,= 0 R./x = 615 X 103/in. 
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TURBULENT BOUNDARY L 
of this paper, while 6; is computed with the modified 
Eq. (37a), 

5, = 0.381 (R,/x)~'*x' 
For R,/x, values of 575,000 per in. and 655,000 per in. 
for M, 1.7 and 2.0, respectively, have been read out of 
It can be seen that the ex- 
For M, 1.7, 
the agreement with the theoretical curve is good as far 


the data of reference 15. 
perimental 5* values follow the x law. 
as the few test points permit a conclusion, while for 
M, 2.0 the test points are about 5 per cent below. 
Correspondingly, the values of momentum thickness 
and friction coefficient determined from the experi- 
mental data are also lower for this Mach Number. 
These statements still are somewhat preliminary be- 
cause, first, the correction for the laminar part involves 
some uncertainty and, second, the data have not been 
corrected for a possible error of pitot tube position. 
Asa first check, however, the agreement is satisfactory. 
For comparison, 6,;* values have been inserted as the 
dotted line in Fig. 11, which are computed with an 
average R,/x = 625,000 per in., and, within +1.5 
per cent, represent the conditions for A/, 1.7 and 2.0 
if compressibility effects are disregarded. 


(VIII) SupPLEMENT 


Since the results of this paper are supported par- 
tially by experiments, the remark appears to be justi- 
fied that similar results can be obtained without resort- 
ing to pipe tests simply by basing the equations for 
incompressible flow [Eq. (37a) and (43a)] on some 
mean values for » and p. There are, of course, several 
possibilities that lead to appreciable differences for 
higher Mach Numbers. The most suitable one will 
have to be decided by experiments. The simplest one 
is to base » and p on the mean static temperature T, 
which can be obtained from Eqs. (4), (7), and (7a) as 


~ | 7 + K(1 — Y’")|dY (55) 
This results in 
T ao 2 _ y= 4 - mn 
T, + pre A =1+4+ “org M,” (56) 
and, form = 7, 
T/Ty = 1 + 0.0444 M,? (56a) 
with 
o(T) py iz /7.3"* 
and 
ul) /ue = (T/Ty)* 


we obtain, for the compressibility effects, from Eqs. 
(37a) and (43a) 


6/6; = (1 + 0.0444 M,2)° +? 


(57) 


and 
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C,/Cr = (1 + 0.0444 M,2)~4-9” = (58) 
At AL, = 5, the ratio 6/6; as obtained from Eq. (57) is 


9 per cent below that of Eq. (38), while the ratio C,/C,; 
as obtained from Eq. (58) at the same Mach Number is 
5 per cent above that of Eq. (44). 


(IX) SUMMARY 


For the turbulent boundary layer on a flat plate in 
compressible flow, the local density is computed with the 
assumptions of constant energy per unit mass, constant 
exponent for the velocity profile, and constant pressure. 
By integration of the local density from the wall to the 
free stream, a term is obtained which is called ‘‘mean 
static density.”’ It is shown that this term is useful for 
rapid determination of different mean values for the 
density from mass and momentum flow at the plate 
and in the pipe, as well as for the ratios of displacement 
thickness and momentum thickness to the full thickness 
of the boundary layer. Numerical results are given in 
graphical form for Mach Numbers 0 to 10 and for ex- 
ponents of the velocity distribution law 5 to 10. 

To obtain the compressibility effect on the boundary- 
layer thickness itself, experimental findings of the fric- 
tion coefficient in pipes are quoted which suggest the 
validity of the Blasius law for compressible flow with 
mean values for density and viscosity based on the as- 
sumption of a mean state. The determination of the 
ratio of mean state temperature to axis temperature 
for pipe flow in connection with this law permits the es- 
tablishment of a relation between the shear stress at 
the wall and conditions at the pipe center. In analogy 
to the principle used first by Prandtl and von Karman 
for incompressible flow, it is assumed that this relation 
can be applied to plate flow, with the free-stream con- 
ditions replacing those at the pipe center and the bound- 
ary thickness replacing the pipe radius. By use of the 


momentum theorem, relations for the variation of 


boundary-layer thickness and friction coefficients with 
Mach Number The results given in 
graphical form are confined to the validity of the '/; 
’* law. 


are obtained. 
power law for velocity distribution and the R, 
With increasing Mach Number, the boundary-layer 
thickness increases slightly and, at 1/, = 10, has about 
The local and mean 
Mach 
at \7, = 10 have about 35 per cent of the incompressible 


twice the incompressible value. 
friction coefficients decrease with Number and 
value. Some estimations are made regarding the prob- 
able range of applicability. Comparison of results 
with those of other theories has evinced considerable 
differences. A subsequent checkup with preliminary 
data from wind-tunnel experiments, however, shows 
satisfactory agreement and supports the validity of the 
assumptions up toa Mach Number of about 2. 
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The Hot-Wire Anemometer in Supersonic Flow 


(Continued from page 572) 


explicit results for velocity and density fluctuations 
and the correlation between the two. 

(5) By using the reasoning given above, it is pos- 
sible to check the validity of a priori assumptions for the 
correlation between various parameters (density-tem- 
perature, velocity-temperature) and obtain further 
check on the validity of the assumption neglecting the 
pressure fluctuation. 
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Turbulent Boundary-Laver Characteristics at 
| Supersonic Speeds—Theory and Experiment’ 


ROBERT 


EK. WILSONt 


The U'nwersity of Texas 


ABSTRACT 


Turbulent boundary-layer characteristics at supersonic speeds 
have been experimentally determined from Pitot pressure surveys 
made at several stations along a flat plate, and the data are com- 
pared with theoretical results. 

The theory is based on the 
differential equation for the velocity distribution in the incom- 
pplies for compressible flow when the 


assumption that von Karman’s 


pressible boundary layer a 
variation in density through the boundary layer is taken into 
Taking the case of'a thermally insulated flat plate and 
a relation between mean turbulent 


account 
neglecting radiation effects, 
skin-friction coefficient and 
ppears to be valid throughout the 


and 


of 


Number is derived, 
test range 


Reynolds 
the relation a 
Mach Numbers and Reynolds Numbers. 

The check between theory and experiment is limited by the 
maximum test values of Mach Number and Reynolds Number, 
which are approximately 2.2 and 19 X 10%, respectively. The 
theory is expected to be valid at Reynolds Numbers higher than 
19 X 10° for Mach Numbers within the test range. However, 
the validity of the theory for Mach Numbers greatly in excess 
of 2.2 should be checked by further experiments. 


SYMBOLS 


y = distance along flat plate, measured from leading edge 


distance normal to plate, measured from surface 


velocity in x direction 

p = density 

T = temperature 

“« = absolute viscosity 

= Mach Number 

R = Reynolds Number 

C; = local skin friction coefficient 

Cy = mean skin friction coefficient 

y = ratioof specific heats = 1.4 for air 
Subscripts 


= conditions at the surface of the plate 


free-stream conditions 


INTRODUCTION 


(I) 
pare OF BOUNDARY-LAYER TESTS have been 
conducted in the Mach Number 1.73, 2.00, and 
2.25 nozzles of the Ordnance Aerophysics Laboratory 
wind tunnel.'!| Measurements were made in the bound- 
ary layer on a flat plate that spanned the closed test 
section of the tunnel. The temperature difference be- 
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tween the plate and the tunnel walls was small so that 
the boundary-layer flow corresponded closely to that 
over a thermally insulated plate and, also, radiation 
effects were negligible. Turbulent velocity profiles and 
mean skin friction coefficients were determined from 
Pitot pressure surveys of the boundary layer made at 
several stations along the:plate. A theory for turbulent 
flow? which takes into account the effect of compressi- 
bility is presented and compared with the experimental 
results. 

RESULTS 


(II) THEORETICAL 


(A) Boundary-Layer Velocity Distribution 


The theoretical treatment is based on the assumption 
that von Karman's differential equation for the velocity 
distribution in the boundary layer 
applies for compressible flow when the variation in 
density through the boundary layer is taken into ac 
Taking von Karman’s equation, 


incompressible 


count. 


Vr/p = —x[(dU/dy)?/(d?U/dy*)| (1) 


the shear stress, 7, is assumed to be constant through 
the boundary layer and equal to the value at the sur 
face of the plate. The von Karman’s ‘“‘uni- 
constant, «, is assumed to be unaffected by com- 
the value that has 

for incompressible 


value of 
versal’ 
and equal to O04, 

determined 


pressibility 
been 
flow. 

An analysis based on Eq. (1) has been made by Frankl 
and Voishel.‘ However, the relations given by them are 


not in a useful form, and some of the results are based 


experimentally 


on mathematical approximations that make them in- 
applicable except at small Mach Numbers. 


After writing 


Vr/pe = U, (2) 
re) U/U, (3) 
n = px U,V/ px (4) 


(1) becomes 


Eq. 
d log (do/dn) |p 


“ do : —e. se 


Taking the case of a thermally insulated flat plate, 
assuming that the 


(d*p/dn?) 
(d@/dn)* 


neglecting radiation effects, and 
energy per unit mass in the boundary layer is constant, 


the density variation is given by 
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Pi 1 | Meal i V2 ( ) (6) Eq. (7) may be readily integrated twice to obtain the ; 


vr. 4 ” r « r.lawe re] +4 r 1c ¢ : > ® 

p . 2 2 U; boundary-layer velocity distribution. Following 
Frankl and Voishel,* the constants of integration ; 

The density at the surface of the plate is obtained , . ieee 
=." ' poate, Po, 0 . evaluated at the edge of the laminar sublayer, where 





by putting U = Oin Eq. (6); and Eqs. (3), (5), and (6) 
are combined to yield n= o=5 
(d@/dn),-; =f 
d log (do@/dn) = 


= (7) The constants s and f are assumed to be unaffected }y 


d V1 — a( 2 sioaes SiGe 
? P/O compressibility and equal to the values for incompres. | 
where sible flow.’ The first integration of Eq. (7) results jy | 
[(y — 1)/2] A? dn/dd = (1 foe” 4 
¢ = . (8) 
1+ [(y — 1)/2] 1° where 


B = (xdi/Ve) {sin—' [Vo (6/d1)| — sin—! [Vo(s/¢,)] } 
Integration of Eq. (9) yields the velocity distribution 


ets oilkdi V 1 — o(¢/ o1)" + o(G/¢)) | oo oi[xgi V 1 — o(s gi)? + a(s/q)| 
fl (ko1)? + o| f[ (xg)? + a] 


Unless the Reynolds Number is extremely low, it can be demonstrated that s < n, except in a small region near 
the wall; and, also, the second term on the right-hand side of Eq. (11) is much smaller than the first term. There 
fore, for the greater portion of the boundary layer, the velocity distribution may be taken as 


_ dilkdiV 1 — o(G/ di)* + o(G/i)1 Os 
fl (ko)? + o| 


When comparing the theoretical velocity distribution with experimental results, it is convenient to plot y 4 ys 
U’/U,, where @ is the boundary-layer momentum thickness. In order to put Eq. (12) in these terms, it is first 





necessary to compute 9, which is defined by the following equation : 


§ = / (pl pl) [1 — (U/U))|dy (13 
“ ( 
where 6 equals boundary-layer thickness. With the use of Eqs. (3), (4), (6), and (8), Eq. (13) may be put in the 


following form: 


px U9 _ 91 he go) (1 — 2/6) (a 


- , |[l — oa (b/¢,)?| do 


1 ( (14) 
bx 1+ [(y — 1)/2] 41,’ , ) _— i 





The derivative dn d@ under the integral in Eq. (14) is given by Eq. (9) for the portion of the boundary layer outside 
the laminar sublayer. However, at high Reynolds Numbers, a negligible error is made in evaluating 4 if Eq. (! 
is assumed to hold for the entire boundary layer. After substituting the value of dy/d@ from Eq. (9) into Eq 
(14), it appears that the integral cannot be evaluated in closed form except for .1/; = 0. In order to obtain at} 
approximate expression for 9, the integral was evaluated numerically over a range of Mach Numbers from 0 to I! 
and over a wide range of values for @. Using the integrated expression for @ at .1/; = 0 as a guide to the form?| 
should take for 1/; > 0, it was found that the following approximate expression checked all the numerical results] 


within a maximum error of about 5 per cent. 


ry ——* 
px LU’, sii o1(Kd 2) elxoi/Ve)tsin-'y/e — sin-'[ Vo(s/d0)] (13 

Mx f(«d1)?* ; 1 ob [(y — 1) 2] M,"} 
By dividing Eq. (12) by Eq. (15), the velocity distribution can be obtained in terms of y/@ and UU). Making 


use of Eqs. (3) and (4), the final result is 


y {E+ [Cy = 2/2) AN} (a6)? Led V1 = 0 U/ 0)? + (U/C) soaorVertsin-AVeit UN) = sin-Val (46 


d (xo: — 2) [(xo,)* + a| 
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An expression for # has already been derived [see Eq. 
(15)|. However, before computing d@ dx, an approxi- 
mation is made which considerably simplifies the re- 
maining calculations. The term sin~! [+/o(s $;)] in 
the exponent of e in Eq. (15) is taken as 


When comparing Eq. (16) with experimental re- 
sults, it is desirable to express ¢; in terms of the local 
sults, . : 
skin friction coefficient, C;,, based on density at the 


surface of the plate. By definition, 


C;, = 7/(1/2) ps. U1? (17) 
“3 1 , . ~ f » 

es sin a(s ;=> a(s ) (20) 

and, from Eqs. (2), (3), and (17), I | ’ m 
: The accuracy of this approximation increases as the 

fi = V2/C (18) 1 , 

= relative thickness of the laminar sublayer decreases, 
(B) Skin Friction Coefficient ind this occur with incre ising Reynolds Number 
By introducing the Reynolds Number based on density 
‘ § ; ) 


By equating the skin friction drag on one surface of and viscosity at the wall, 
the plate to the momentum deficit in the boundary 
laver and then differentiating, it can be shown that the 
hear stress at the wall and the momentum thickness , saa , z = 
sheer sires and by making use of Eqs. (6), (15), (17), (1S), and 


Rx = px Ux Mx (21) 


are related by (20), it can be shown that Eq. (19) may be written 
tT = p,U,?(d0/dx) (19) in the following form: 
2 sin! 2 sin !~/a 2x" 
fre c “fo Ff C;’ Va Cy, 


After integrating Eq. (22) and evaluating the constant of integration for the case of turbulent flow over the 
co when Ry = 0, the Reynolds Number is given in terms of the local skin 


entire surface of the plate so that C;, = 


friction coefficient as follows: 


2 [x Vo 2 Vo = v2 
* fixie” é ( sin~'/o \ ry =k sin 'VYo sin~!+/o@ 


9) j 

wh (2 + fs ) (23) 
. KS 1 f/ — e 1 -* 
jee’ sn "vo sin 'WYo 


Following von Karman’s derivation for the case of incompressible flow at high Reynolds Numbers,’ all the terms 
on the right-hand side of Eq. (23) are neglected except the one containing 1/C;, so that the relation between local 


skin friction coeflicient and Reynolds Number is given by 


in ! fre" 
(= \ *) == log (” ) + log (Cy, Re) (24) 
Vo VC, KV/2 4 kvV/2 


The mean turbulent skin friction coeflicient based on 


It should be noted that some of the terms in Eq. (23) 
density at the surface of the plate may be obtained 


which were discarded, although they are smaller than 


the term kept, are not negligible. This fact is taken from 

into account in the case of incompressible flow by ad- D _ 

justing the constants of the theoretical relation to con- Cr, = — = Cy,dRx (26) 
(i/Z ) Px [ 4°X Ry 0 


form with experimental data. 

With the assumption that the constants «, s, and f 
are unaffected by compressibility, the first term on the 
right-hand side of Eq. (24) and the coefficient of 
log (C;,R,) may be taken from the incompressible flow  cvefficient obtained from Eq. (26) for incompressible 
flow can, by making certain approximations, be written 
in a form identical to that for local skin friction co- 
efficient with somewhat different values for the nu- 
The same procedure holds here, 


where D equals the skin friction drag on one side of the 
plate from the leading edge to station x. According 
to Schoenherr,® the expression for mean skin friction 


case. Taking von Karman’s numerical values,’ Eq. 


(24) becomes 


(= la/o 
Vo 


At M, = 0, the above relation reduces to von Karman’s 


] ‘cal constants 

= a . nad merical constants. 

f— = 1 + 4.15 logiy (Cy, Re) (25) : : J 

i with the numerical constants again unchanged by com- 
pressibility. Taking Schoenherr’s values for the con- 
stants,® the relation between mean skin friction coeffi- 


result for incompressible flow. cient and Reynolds Number is given by 
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(sin-!+/¢/+/oe) (0.242 V Cr,) = logw (Cr, Re) (27) 


Relations between skin friction coefficient and Rey- 
nolds Number based on free-stream conditions may be 
obtained from Eqs. (25) and (27). Brainerd and 
Emmons’ have shown that an approximate relation 
between viscosity and temperature is 


bx /pi = (T%/T))9-7% 


Taking the above expression, it can easily be shown 
that 


pill 1X 


a _ ae R. 
Mi T; 

: us 7*) Pa 
Cc, = = Cc 

' (1/2)p.U;" & ™ 


D 7) ae 
Cy = = Cs 
5 (1 2) p, Ux (= 


So far, it has been assumed that the flow in the bound- 
ary layer is isoenergetic, and this assumption results 
in the relation, obtained from Eq. (6), 


R= 


(28) 


T/T, = 1+ [(y — 1)/2]M,? (29) 
The assumption is in error, even for a thermally in- 
sulated plate, if the Prandtl Number is different from 
unity, as in the case of air; and, before proceeding with 
the calculations, the error is briefly discussed. Taking 
the case of a thermally insulated plate and neglecting 
the effects of radiation, calculations have been made! 
which indicate that the assumption of isoenergetic 
flow should have a negligible effect on the accuracy of 
the calculation for momentum thickness. This re- 
sult is due to the fact that, while the flow is not iso- 
energetic, for the case of the thermally insulated plate 
there is no net energy deficit or increase in the boundary 
layer. Thus it appears that the assumption of iso- 
energetic flow should have a negligible effect on the 
factor sin~!4/o/+/o in the skin friction relations, since 
this factor resulted from integrating across the bound- 
ary layer to determine the momentum thickness. 
But it does seem reasonable to take the error in the 
assumption into account by basing the skin friction 
coefficients and Reynolds Numbers in Eqs. (25) and 
(27) on the actual values of density and viscosity at the 
wall. Therefore, rather than use Eq. (29), the tem- 
perature ratio is taken as 

T/T, = 1+ G[(y — 1)/2] M7 (30) 
where C, equals Hilton's temperature recovery co- 
efficient.* For subsonic speeds, when the effect of 
compressibility is not large, C, has been theoretically 
computed by Squire* as the cube root of Prandtl’s 
Number. From Eqs. (25), (27), (28), and (30), the 
skin friction coefficient relations based on free-stream 


conditions are given by 
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sin ~!+/o¢ _y-!1 : ae 
= ( l + C 6 M, ~ o = 
Vo 2 Wf, 


LZ + 4.15 logo (C,R) —_ 


— | 
3.19 logw (1 +. C, Ls : ue) (31 


—'72 0.242 
M,? ee a 
V Cy 


sin 


l — 
Ye (1467 
Vo 2 


=— J 
logio(CrR) —_ 0.768 logio (: oa Ge L 3 ur) (32 


(C) Comparison with von Karman’s Result 


von Karman' has given an approximate relation | 


for mean turbulent skin friction coefficient as fol 
lows: 
( y—1 - ) 2 0.242 ;, » 
l + / = 10210 (Cr — 
5 1 V Cp $1 KF 





l y¥-—1 
5 logio (1 + = MM] 


The above equation is the result of assuming that the 
flow in the boundary layer is isoenergetic and that the f 
compressible flow relation is identical to that for in-} 
compressible flow if the density and viscosity at the} 
surface of the plate are used in computing skin frictior 
coefficient and Reynolds Number. The factor 1? 
rather than 0.768 appears in the second term on the 
right-hand side of the above equation because of the] 
relation between temperature and viscosity used by vot} 
Karman. Thus it is seen that the most significant con | 
tribution of the calculations given here is the introduc 
tion of the factor sin-'V/o/~/o. Taking an exper 
mental value of 0.88 for C,,* Eq. (32) is compared with} 
von Karman’s result for a wide range of Reynolds 
Numbers and Mach Numbers on Fig. 1. The preset! 
theory gives skin friction coefficients that are consider } 
ably higher than those given by von Karman’s result} 
principally because of the factor sin~!'Y/oa/ Vo. ; 
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(III) COMPARISON OF THEORY AND EXPERIMENT 


(A) Velocity Distribution 

Pitot pressure surveys have been made of the bound- 
ary 
The tests were run in the Ordnance Aerophysics Lab- 
oratory wind tunnel using the Mach Number 1.73, 
2.00, and 2.25 nozzles, with some additional variation 
in the Mach Number of the flow over the plate being 
obtained by changing angle of attack. A sketch of 
the plate used for the experiments is shown on Fig. 2, 
and the Pitot tube is shown on Fig. 3. In addition to 
the flat plate tests, measurements were made in the 
boundary layer on the tunnel wall to determine the 
effect of aerodynamic interference of the Pitot tube on 
the experimental velocity profiles."! 

An experimental velocity profile from the measure- 
ments made on the ttinnel wall has been chosen for 
comparison on Fig. 4 with the theoretical result given 
by Eq. (16). This profile was chosen for the compari- 
son rather than a profile measured on the flat plate 
because, using the same diameter tube, Pitot pressures 
could be measured relatively closer to the surface in the 
tunnel wall tests. The value of Cy, given on Fig. 4 
was not measured but was obtained indirectly from the 
measured momentum thickness of the boundary layer 
and the theoretical relations for skin friction coeffi- 
cient. There are appreciable differences between the 
theory and experiment, and the differences noted on 
Fig. 4 are essentially the same for velocity profiles meas- 
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ured on the flat plate. This result is not surprising, 
since, if the velocity profile for the incompressible flow 
case derived from Eq. (1) is compared with the experi- 
mental results of Zijnen,'’ it is found that there are 
disagreements of the same order of magnitude as those 
shown on Fig. 4. 

Prior to the derivation of Eq. (1), both von Karman 
and Prandtl made calculations for incompressible flow 
using an expression for the velocity distribution of the 


form 


U/U, = (y/6)'™ (33) 
Taking a value of 7 for n, Zijnen obtained good agree- 
ment between Eq. (33) and experimental results for in 
compressible flow over a moderate range of Reynolds 
Numbers.'? However, Prandtl'* has indicated that 
n is a variable and increases with Reynolds Number. 
The experimental results for compressible flow are simi- 
lar to those for the incompressible flow case. It has 
been shown that, using m = 7, Eq. (335) agrees well with 
the experimental data in the test range of Mach Num- 
bers over the flat plate from 1.6 to 2.2 for Reynolds 
Numbers from about 3 xX 10° to 19 XK 10°! For 
higher Reynolds Numbers, » should be increased, and 
on Fig. 4 it is shown that excellent agreement with the 
velocity distribution on the tunnel wall is obtained 
using m = 9. The “equivalent” flat plate Reynolds 
Number of the profile shown on Fig. + was computed 
from the measured momentum thickness and was found 
to be 50 X 10°. 


(B) Shape Parameter 


When computing the boundary-layer shape param- 
eter, //, defined as the ratio of displacement thickness 
to momentum thickness, Eq. (33) should be used be- 
cause of the agreement with experiment. If 7 is taken 
as 7, the assumption of isoenergetic flow is used in com- 
puting momentum thickness from Eq. (13), and the 
displacement thickness is computed from 


bw) 
| [1 — (pU 
0 


then the shape parameter is given by 


Sk 


6* = pil)|dy 


1+ {(y — 10/2) 40? o! 
a 7 (8 + 28) 
a x go’ 


u (S + 22) (9 + 22) 


(34) 


Values of 7 calculated from the above expression are 
compared with the experimental results on Fig. 5, and 
the agreement As already mentioned, 
the value of m in Eq. (33) varies somewhat with Rey- 
nolds Number, and this fact accounts, to a consider- 
able extent, for the slight scatter of the data on Fig. 
». Fortunately the value of // is not very sensitive to 


is excellent. 
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FIG. 5 SHAPE PARAMETER VS. MACH NUMBER 
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n, and it is believed that Eq. (34) should be satisfactory 
for most purposes. 
(C) Skin Friction Coefficient 

Experimental values of mean skin friction coefficient 
were determined from the results of the Pitot pressure 
surveys on the flat plate by means of the momentum 
equation and the assumption of isoenergetic flow in the 
boundary layer. The skin 
coefficient and Reynolds Number were complicated bj 
the fact that some variation in the Mach Number of 
the flow over the plate was always present. Both 
skin friction coefficient and Reynolds Number were 
based on the average Mach Number, and the details of 
the data analysis have been fully discussed in previous 


calculations of friction 


» 


reports." * 

The test values of mean skin friction coefficient and 
Reynolds Number are plotted on Fig. 6, and it is appar- 
ent that in some cases the flow was laminar for a con- 
siderable distance from the plate leading edge. Using 
the laminar flow results of Brainerd and Emmons,’ 
it is seen that the transition Reynolds Number ex 
ceeded 4 X 10° for the two highest Mach Numbers. 
In order to compare theory with experiment, it 1s 
necessary to correct the test values of Cy and R to 
100 per cent turbulent flow. When the measurements 
are made in a region where turbulent flow is fully de- 
veloped, the frictional force acting on the portion of the 
plate from the leading edge to a particular measuring 
station is taken equal to the frictional force on a plate 
with 100 per cent turbulent flow and with an ‘‘effective 
leading edge located at a distance x)» downstream from 
the actual plate leading edge. For the following calcu- 
lations, the uncorrected test values are denoted by 


Cy’ and R.,. 








Cy and R for 100 per cent turbulent flow 
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are given in terms of Cp’ and R, by 


Cr = [R./(Rz — Rx] Cr’ 


Cai ~ 3 =) 


where R,, equals Reynolds Number based on the dis- 
tance x». After substituting Eqs. (35) into Eq. (32), 
the following theoretical relation between Cy’ and R, 
is obtained: 


R aoe R f M;,) " 
Zz Zo = 1 y 9 ¢ i) ae 
Ye vo’ ~ oe 


— | 
0.768 logi« (1 + C, z o un) (36) 


where 


Sj l if _ 1 ow ise 
f(Mh) = 0.242 (=. = "\( +c 7 us) 
o 4 


(37) 


Eq. (36) is used to compare the theory with test results 
in the region where turbulent flow is fully developed. 
The most direct method to make the comparison would 
be to use measured values of R,, in Eq. (36). The ef- 
fective leading edge is at a station somewhat ahead of 
the transition point so that values of R,, could perhaps 
be deduced from measurements of transition Reynolds 
Number. Transition Reynolds Numbers were not 
measured, and they cannot be determined with any 
degree of accuracy for most of the test Mach Numbers 
by direct extrapolation of the measured skin friction 
coefficients to the theoretical curve for laminar flow. 
It can be seen from the data plotted on Fig. 6 that the 


0.242 
he aac, , et 
 o( > uy) (sin-! + | C; 


logy | 


It may be seen from Eq. (38) that, if Cp and R are multi- 
plied by certain functions, the same curve serves for all 
Mach Numbers. The test data are plotted in this 
fashion on Fig. 8 and are compared with the theoretical 
curves. It should be noted that C, was taken as 0.SS 
in all of the above calculations. 


It is of considerable interest that a good check with 
experimental values of mean skin friction coefficient 
is obtained from the theory presented herein; thus, 
for engineering purposes, it may not be necessary to 
resort to a considerably more complicated theory such 
as that given by Ferrari.'' There appears to be no 
reason why the theoretical relation for mean turbulent 


test Reynolds Numbers were not low enough to perm 
obtaining the transition Reynolds Number by thi 
method except for the two highest Mach Numbers 
Consequently, Eq. (36) has been applied to the tey 
data, and both R,, and f(.\/;) have been determine, 
by the method of least squares. The second term oy 
the right-hand side of Eq. (36) also could have bee 
determined from the data, but the calculations woylg 
have been much more complicated. Since this tery 
is small compared to logyw(Cy’R,), it was computed be 
fore analyzing the data. It is believed that agreemen; 
between values of f(.1/,) obtained from the data in this 
fashion and values computed from Eq. (37) should 
indicate that the theory is reliable. 

The curves given by Eq. (36), with R,, and f(.\/)) de 
termined from the data analysis, are plotted through 
the uncorrected test values of skin friction coefficient 
on Fig. 6. The values of R,, were used in Eqs. (35) t 
correct the data to 100 per cent turbulent flow, and the 
corrected data are compared on Fig. 6 with theoretical 
curves given by Eq. (32). The values of f(.\/;) ob 
tained from the data analysis are compared on Fig. 7 


> 


with theoretical values computed from Eq. (33 
The agreement between theory and experiment on both 
Figs. 6 and 7 is within the scatter of the test data. 

To summarize the results, the test data corrected t 
100 per cent turbulent flow have been plotted in a 
manner such that all the data should fall on the same 
curve regardless of Mach Number. In order to ac 
complish this, it was noted that Eq. (32) can be written 


as follows: 


onl us) =e _R|| 


9 , — i 1.768 (3 
7 é F o( 1 i Cy : ue) | 
l \ a)- ») 


skin friction coefficient should not be valid for Rey 
nolds Numbers outside the test range. However, be 
cause of the various approximations and assumptions 
made in the theory, the validity of the theory at 
Mach Numbers greatly in excess of 2.2 should be 
checked by further experiments. 


(IV) CONCLUSIONS 


(1) Using von Karman's differential equation for 
the boundary-layer velocity distribution in the case 
of incompressible flow, a relation between mean turbu 
lent skin friction coefficient and Reynolds Number for 
compressible flow has been derived and appears to be 
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actory at Mach Numbers up to 2.2 for the case of 


satis! 
, thermally insulated flat plate when radiation effects 


are negligible. 

(2) The reliability of the theoretical values of skin 
friction coefficient for Mach Numbers above 2.2 should 
be checked by further experiments. 

(3) There are discrepancies between the theoretical 
and experimental boundary-layer velocity distributions 
for compressible flow which are of the same order of 
magnitude as the discrepancies between von Karman’s 
theoretical result for two-dimensional incompressible 


flow and experiment. 

(4) The boundary-layer shape parameter, /7, may be 
calculated from a velocity distribution of the form 
U/U, = y 5)'’", where, for most purposes, 7 may be 
assumed to be independent of Reynolds Number and 


Mach Number and equal to 7. 
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RIEF REPORTS of investigations in the aeronautical sciences 

and discussions of papers published in the JOURNAL are 

presented in this spec ial department. The publication will be com- 

pleted approximately 8 weeks after receipt of the material. The 

Editorial Committee does not hold itself responsible for the opinions 

expressed by the correspondents. Contributions should not exceed 
800 words in length 


Comments on ‘‘Heat Transfer in Laminar 
Compressible Boundary Layer. . .”’ 


PD, Grootenhuis 
Mechanical Engineering Department, Imperial College of Science 
and Technology, London, England 


May 26, 1950 


\ THE PAPER ‘“‘Heat Transfer in Laminar Compressible Bound- 
he Layer on a Porous Flat Plate with Fluid Injection” 
JOURNAL OF THE AERONAUTICAL SCIENCES, Vol. 16, No. 12, p. 
741, December, 1949) S. W. Yuan states that the amount of 
heat removed by the plate per unit area is equal to the heat ab- 


sorbed by the coolant, which leads to Eq. (20) 


*L ve 
o7 
/ al ) dx = putolp (Tw — To)L 
0 3 fie 


This equation integrates the heat conducted to the plate through 
the boundary layer over the length Z of the flat plate (left-hand 
side expression) and equates this to the heat absorbed by the 
coolant of inlet temperature 7) on the assumption that the cool 
nt attains the temperature of the plate 7, at all points. It 
should be noticed that the left-hand side of the equation is 
treating the temperature gra- 


integrated over the limits O-L, 
dients in the boundary layer, reckoned for the surface conditions 
as a variable of x, whereas the surface temperature of the plate 
is taken to be independent of x. This cannot be so for the initial 
conditions assumed. The analysis in the paper by Yuan as- 
sumes a constant velocity of injection along the flat plate. This 
leads to a decreasing surface temperature with increase of the 
distance x from the leading edge of the plate. Consequently, 
the basis of Eq. (20) is in error, and the subsequent estimates 
of the quantity of coolants required will have to be revised 

That the surface temperature decreases along the plate for 
a uniform velocity of injection has been shown by the writer in 
his joint paper with N. P. W. Moore ‘“‘Some Observations on the 
Mechanism of Sweat Cooling,’’ published in the Proceedings of 
the 7th International Congress of Applied Mechanics (London, 
1948, Vol. 3, p. 106). 

A further inaccuracy is introduced by assuming that the coolant 
attains the surface temperature as it is forced through the porous 
material. This need not always be the case. The temperature 
rise of the fluid on its passage through the porous material de 
pends on the internal heat-transfer coefficient and the thermal 
conductivity and porosity of this material. The theoretical 
analysis of the process involved and an estimation of the tem 
perature gradients in the porous solid and in the coolant were also 
given in the writer’s paper already. referred to. This showed 
that the temperature of the coolant on emerging at the hot sur- 
face may, under certain conditions, be considerably less than the 
surface temperature 
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Comment on ‘‘Heat Source in a Uniform 
Flow’’* 


B. L. Hicks, W. H. Hebrank, and S. Kravitz 
Ballistic Reszarch Laboratories, Aberdeen Proving Ground, Md 
June 9, 1959 


W'° SHOULD LIKE to point out the relationship between our 
theoretical investigations and those of Tsien and Beilock 
In our early work, we had described the wake downstream of a 
heated region and shown the connection between the vorticity 
gradient along the stream and the gradient of heat addition nor- 
mal to the flow.' In later work,? we examined various types of 
diabatic flow including that in the neighborhood of a continuously 
distributed collection of heat sources proportional toe~®”. Our 
considerations here differed from those of Tsien and Beilock prin 
cipally in the fact that we did not consider supersonic flow or 
pass to the limit of a point heat source but considered one ex 
ample in some detail. It can be seen from our calculations, for 
example, how the source of heat acted like a source of fluid, could 
decrease the fluid velocity, and even produce stagnation points if 
the heat source were strong enough. 

Several defects of the theory in reference 2 have now been re- 
moved.’ Comparison of the results of that paper with the re- 
sults of Tsien and Beilock showed us that there had been a slight 
error in our application of the Glauert-Prandt] approximation for 
compressibility of the fluid. Also, the theory in the symposium 
paper was based on a first-order perturbation calculation that 
assumed that the vorticity was negligible. We have since shown 
that the vorticity and also the momentum transport (p + 
pv”) are indeed of the second order in the perturbation velocities 
and their derivatives. 

Moreover, it was still evident that in reference 2 some error was 
inherent in the perturbation calculations because the stagnation 


pressure passed through a minimum. It is known that a perfect 


* Tsien, H. S., and Beiilock, M., Readers’ Forum, Journal of the Aero- 
nautical Sciences, Vol. 16, No. 12, p. 756, December, 1949 
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gas containing only positive sources of heat cannot exhibit in- 
crease of stagnation pressure. This anomaly has been resolved 
by a more careful analysis of the density perturbation along the 
axis of the flow. The results of this analysis are shown in Fig. 1. 
The values for curve I (and for Fig. 5b of reference 2) were 
plotted from 


Mel(V\2 
Pw ty? —)(14+"%)-1 
Po Po 2 J 0 Po 


Values for curve II were plotted from the results of numerically 
integrating the almost exact equation [see reference 1, Eq. (2.5)] 


a In —2yMo V1 — Me? ; 
Pr oe , ge (1 — W) xX 
Ox Vi7y- 1)/2 
To\? — 1) (1 — eo" 
1+ (y q1 é 
7 Mo Lal 


in which the Crocco number can be calculated from the equa- 
tions derived in reference 2 for the velocity and pressure vari- 
ations and 
qe = [Moe? + 2/(y — 1))7?¢0 

Comparison of curves I and II for go = 0.0088, corresponding 
to about 10 per cent stoichiometric, show a large inaccuracy due 
As a further check, the 
calculations were repeated for go = 0.001, which is small enough so 


to density perturbation along the x axis. 


that the equation for curve I should be reasonably accurate. 
From Fig. 1 it can be seen that, for the small go value, the curves 
I and II do lie close together. 

It should be emphasized that the velocity and pressure per- 
turbations, in the case of heat source in a uniform flow, can be 
reasonably small, as in the above calculations, while the density 
and temperature perturbations are larger than is admissible in a 
We suggest, therefore, that an attempt be 
made to construct a new theory that makes allowance for large 


first-order theory. 


variation of quantities along the flow direction (as calculated, for 
example, from well known one-dimensional diabatic flow theory), 
but still assumes small variation of quantities in directions nor- 
mal to the flow. This would amount to modifying one-dimen- 
sional theory by allowing for some sidewise expansion. 

We should also like to call attention to a calculation of the force 
exerted on a cylinder in a field of diabatic flow.4 Compressibility 
corrections have been made, and it appears that the calculation 
can be extended to the case of two neighboring bodies. 

We should like to express our thanks to Professor Tsien for 


sending us a prepublication copy of his paper. 
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On a Variational Principle in Lifting-Line 
Theory 


Alexander H. Flax 
Cornell Aeronautical Laboratory, Inc., Buffalo, N.Y. 
June 2, 1950 


_— PROBLEM OF DETERMINING the lift and induced drag dis- 
tributions on lifting elements of finite span, based upon 
Prandtl’s lifting-line theory, has been studied by a large number 
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of investigators. 
the development of techniques for the solution of the 
equation of lifting-line theory. 


finite span monoplane in uniform flow, several highly succesc 
4 


Essentially, this problem resolves itself int 
integral 
In particular, for the case of the 
ioe ul 
In addition, there js avai 
able the empirical method of Schrenk,® which provides, wit, 
relatively little computation, startlingly good approximatioy 


methods of solution are available.! 


to the more exact solutions for many cases. In view of these 
facts, it might appear that the presentation of yet another ap 
proach to the solution of so old a problem is hardly worth whi 
Nevertheless, in the opinion of the writer, the variational] prin 
ciple to be discussed here offers a sufficient number of possibj, 
advantages, both in its potentialities for solving certain specifi 
problems and in providing insight into and unifying the genera) 
theory, to make it worthy of consideration. 

The variational principle that has been evolved for lifting-ling 
theory replaces the integral equation of lifting-line theory by th; 
variational equation 


AS (er — doa)? (c/ao) dx + S cai cdx|] = 0 | 
or 
67, + Iz) = 0 9 
where 
C] = the local lift coefficient 


Caz| =c,)(w/V)] = the local induced drag coefficient 


ao = the local slope of the lift curve 

x = the coordinate along the span 

w = the local induced velocity normal to the 
surface 

V = the free-stream velocity 


By a simple application of the calculus of variations,’ it is easily 
seen that this equation is equivalent to the usual integral equa- 
tion of lifting-line theory. 
straightforward, giving 


The variation of the first integral is 


6, = S Aer — aya) (c/ao) dc, dx (3 


The evaluation of the variation of the second integral may be 
readily accomplished by the use of Munk’s stagger theorem,’ 
which states that, if an element of lift is added to a lifting system, 
the induced drag produced is independent of the chordwise posi 


tion of the element. Consider now 


6], = 6S cw V)e dx (4 
If the element of lift involved in the variational process is added 


far downstream of the wing, in the Trefftz plane, the induced 
velocity will be twice that at the wing, so that 


bh, = 2S (w V)éc;c dx 5 
Thus, Eq. (1) reduces to 
J 2[ (cr — ava) (1/ao) + (w/V) icy dx = 0 6 


Since 6c; is arbitrary, this gives 
] 


cj = adola -— te V)] ‘ 


which is the usual lifting-line equation, when for w there is sub- 
stituted 


w= (1/4r) J I'() [dé/(x — £8)] 8 
and 
r= ¢cvV/2 (9 


The same result can be obtained straightforwardly, without the 
use of the Munk stagger theorem, by using for J2 the form of the 
induced drag integral for straight monoplane wings, 


I= 2f fv'(er'(x) log |x — & 


dé dx (10 
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READERS’ 


vhich can be obtained from Eqs. (4), (8), and (9) by integration 
, 


by parts. Ae __ : 
From Eq. (1), it may be seen that the variational principle is 


actually also a minimum principle, since the first integral, having 
a squared integrand, is obviously positive definite, while the 
second integral, representing induced drag, is also positive defi 

nite. Thus, the principle states that the correct lift distribution 
is the one that minimizes the induced drag plus the first integral, 
which can be interpreted as the induced drag of an ideal auxiliary 
wing of the same span as the actual wing, itself giving rise to no 
induced drag, loaded at every station with a lift, Ac,, equal to the 
decrease in lift of the actual wing from the two-dimensional value, 


and working in the downwash field at the actual wing. This 
interpretation follows from the relations 

Ac; = dha — C) (11) 

w/V = (aoa — C;)/ao (12) 

(1:3) 


i= r i Ac)(w/V)e dx 


The practical value of a minimum principle of this type is that 
it lends itself well to approximations, particularly to the Ray 
leigh-Ritz method. In this case, the approach is to assume the 
lift distribution to be given by a series 


C= Yan fa(x) 


where the f,(x)’s are suitable arbitrary functions satisfying the 
boundary conditions on c,; at the tips and leading to finite values 
of the induced drag integrals. The variational condition (2) is 


then replaced by 


(14) 


(15) 


(0/da,) (1, + Is) = O 


The number of terms, m, chosen is then dependent on the accu- 
racy desired. In the theory of the monoplane wing, if the f,(x)’s 
are taken as trigonometric functions sin ”@, where x = —s cos 6 
as in the Glauert! and Lotz? methods, Eq. (15) leads directly 
to equations analogous to the Lotz method. In fact, they be- 
come identical to the Lotz equations if the same series representa 
tions for the chord and angle distributions are also used. If the 
Sears‘ eigenfunctions are employed for the f,,(x)’s, then Eq. (15) 
gives the eigenfunction obtained by Sears. The 
Schrenk method,® from the point of view of Eq. (1), consists of 
averaging a function that minimizes J, (actually this integral is 
then zero) with one that minimizes J., subject to the condition of 
Thus, a number of seemingly unrelated methods 


expansion 


a given lift. 
can be seen to be related in terms of Eq. (1). 

It has been found that the variational principle can be extended 
to more general cases, many of which have not been satisfactorily 
solved. In particular, the cases of wings on a cylindrical body, 
wings with finite end plates, wings with bends in a plane normal 
to the flow, and wings in a field of nonuniform spanwise velocity 
can be treated. In many of these cases, the lift distribution for 
minimum induced drag is already known, thus providing a good 
arbitrary function, fi(x), for the series corresponding to Eq. (14) 
By choosing other functions that minimize J, alone (without, 
however, making J, infinite), it appears that good approximate 
results could be obtained with relatively few terms by the use of 
the variational principle. A forthcoming paper will present a 
more detailed exposition of the variational method in lifting-line 
theory, particularly with regard to these more general cases. 
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Lift of Planing Surfaces* 


B. V. Korvin-Kroukovsky 
Research Professor in Fluid Dynamics, Experimental Towing Tank, 
Stevens Institute of Technology, Hoboken, N.J 


June 15, 1950 


6 bow NOTE IS PREPARED to show that the lift of a flat plate 
planing on the surface of water at an angle of trim r, as pre 
dicted by the mathematical analysis, agrees with the experi 
mental data if proper account is taken of the downwash velocity, 
which is induced by the lift force in the case of a planing surface 
of finite beam. 

In reference 1, the analytical solution for the potential flow 
about a planing surface, originally developed by Herbert Wagner, 
is presented in detail by Pierson and Leshnover. An inviscid 
weightless fluid is assumed, and the planing surface is taken to 
be infinite in span and in its forward length, along which the 
spray sheet flows, and to have a finite trailing edge from which 
the bulk of the fluid breaks away smoothly. The normal force 
acting on unit width of the infinite span is shown to be 


. l , 27 sin a 
F = pli ? 
4 l — COS a . 
1 + cos a — (1 — cos a) In + mrsina 
2 cos a 
(1) 

where 

p = mass density 

V = velocity 

a = angle of incidence 

/ = wetted length, defined as the distance from the trailing 


edge to the perpendicular to the planing surface, drawn 
tangent to the free surface of the spray root 


The discussion will be limited to sufficiently small angles a so 
that cos a can be taken as | and the total lift Z, acting on the 
width 4 of the planing surface of infinite span, can be taken as 
For these small values of a, 


equal to the normal force OF. 
\(1 — cos a)/2 cos] is fractional so that the logarithm is negative 


and the whole term 


—(1 — cos a) In [(1 — cos a)/2 cos a 
is positive but very small so that it may be neglected. To this 
approximation then 
(l + cosa) = 2 
and Eq. (1), rewritten in coefficient forin is 
. g, Qn sin a 
Cc; = =f(a) = 2) 


(1 2)pbl V? 2+ sin a 

The lift force can be considered as the reaction to the vertical 
momentum imparted to the water per second—i.e., L = mw, 
where ™ is the time rate of change of the virtual mass, taken ‘in 
this case to be the mass of a half-cylinder of water of diameter 
equal to the beam of the surface, 6, and of length, V, and where 
w is the vertical velocity imparted to the water. The lift is ex- 


pressed, therefore, as 


* The material reported above was obtained in the course of the Re 
Planing Surfaces at the Experimental Towing Tank of the 


search on 
of Technology under Contract with the Office of Naval 


Stevens Institute 


Research 
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L = (prb?/8)Vw 


Consideration of the energy imparted to the water indicates 
that at the center of pressure of the planing surface half the ver 
tical velocity w is imparted to the water and the angle of down- 
wash is found to be 


e = w/2V = 4L/prb?V? (3) 


as given by Wagner.4 Expressing Z in terms of C; on the basis of 
Eq. (2), 


21 2r 
e= Ci = Ci (4) 
rb 7 
or, expressing e in degrees, 
e = 36.5AC) (5) 


Eq. (2) can be considered valid for a planing surface of finite 
span b, provided the angle a is measured between the planing sur- 
face and the velocity vector at its center of pressure, which is in- 
clined by the angle ¢ to the undisturbed water. The angle 7 be- 
tween the planing surface and the undisturbed water surface is a 
sum of two angles, 


ate (6) 


Table 1 gives values of Ci, e, and 7, for several values of a. 
In comparison with the airfoil data, the values of « are extremely 
large, partly because the half-cylinder is used in computing the 
virtual mass for the planing surface case instead of the complete 
cylinder used for airfoils or deeply submerged hydrofoils and 
partly because of the large values of the length/beam ratio 
(small aspect ratio) used for planing surfaces. 


T = 


TABLE | 
Computation of the Downwash Velocity and Angle of Trim for 


Flat Planing Surfaces 

’ -e« (Deg.) . r= a +e (Dez.)— 
Deg 7 sina Cl r 1 A=2 A=3 A=1aA=2 A 3 
0.25 0.0137 0.0136 0.50 0.99 1.49 0.75 1.24 1.74 
0.50 0.0274 0.0270 0.98 1.97 2.95 148 2.47 3.45 
1.0 0.0550 0.0535 1.95 3.91 5.86 2.95 +.91 6.86 
1.5 0.0822 0.0792 2.89 5.78 8.67 4.39 7.28 10.17 
2.0 0.1095 0.1088 3.79 7.68 11.40 5.79 9.68 13.40 
2.5 0.1370 0.1282 4.68 9.36 14.04 7.18 11.86 16.54 
3.0 0.1642 0.1518 11 O08 16.65 8.54 14.08 19.65 
3.5 0.1915 0.175 12.76 9.88 16.26 

4 0.219 0.1975 14.42 11.21 18.42 

5 0.274 0.241 13.8 

6 0.328 0.282 16.3 

7 0.383 0.322 183.75 





The fact that @ is small compared to 7 permits further simpli- 
fication of Eq. (2). Neglecting z sin a in comparison with 2 in 
the denominator, Eq. (2) may be written as 

C; = 


rsina = 


Ta 
and 
dCi/da = xr (7) 
Differentiating Eq. (6) with respect to C; and substituting 
da/dC; = 1 2\/m [from Eq. (4)] leads to the 


expression 


nm and de/dC; = 


dC./dr = w/(1 + 2d) 


which compares with the form 


4+- (1/A)] 


dCi/da = xw/{(1/2) 


given in airfoil theory, since \ = 1/4, and since half of the cylin- 


der is considered in the virtual mass of a planing surface. 


For 7 in degrees, the numerator of Eq. (8) becomes 7/57.3 = 
0.055 
the ultimate limiting value corresponding to the most favorable 
lift distribution along the beam and, in practice, is found to be 


However, the value of dC;/dr given by this coefficient is 


Thus, for airfoils, as used on actual airplanes, it is cus- 
0.110. It 


smaller. 


to:nary to take 0.10 in place of 27/57.3 = appears 
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that for planing surfaces the best correlation with €xperimenty) 
data is obtained by using the figure 0.04 and writing 


dCi/dr = 0.04/(1 + 2d) per deg 





Eqs. (8) and (9) apply to a surface of large aspect ratio, th 
action of which is approximated by a lifting line Physically, j 
takes into account the vertical deflection of the stream of Water 
by the action of the lift force but does not consider the laterg 
flows present at a surface of a small span.  Siler® demonstrat, 
that, in the case of a surface of small aspect ratio, the effects 
longitudinal and of transverse flows can be simply added. fy, 
a surface of infinite length, finite beam +, in a stream of velocity 
V and angle of trim 7, the normal force per unit length dye 
the transverse component of the flow, V sin r, is given by Lamb x 


P = pbV?[xr/(4 + w)] sin? + 


In coefficient form, denoting by AC, the contribution to ¢,}y 
the transverse flow, 


Pl 2r 





AC; = ~ >, = sin? 7 > 0.887? 
(1/2)pbl V? 4+7 
and 
dACi/dr = 1.76r 
or, for 7 in degrees, i 
; = | 
dACi/dr = 0.000547 1] 
Now, combining Eqs. (9) and (11) gives 
dCi/dr = [0.04/(1 + 2d)] + 0.000547 12 
The values of dC;/dr computed from Eq. (12) for the angk 
7r = 5° are given on the first line in Table 2. On the third line, 
the values of C; for 7 = 10° are given, as obtained by integration 
of Eq. (12), 
Cr = [0.04/(1 + 2d)]7r + 0.000277? (13 
TABLE 2 
Comparison of Analytical and Empirical Values of dC;/dr and ( 
r 1/b 
0.5 1 2 3 1 5 | 
dCi/dr for r 5 ' 
(1) Analytical Eq. (12) 0.0226 0.0159 0.0106 0.0083 0.0070 0.006 i 
(2) Empirical Eq. (15 0.0220 0.0155 0.0110 0.0090 0.0078 0.0070 | 
Ciforr = 10 ' 
(3) Analytical Eq. (13) 0.226 0.159 0.106 0.083 0.070 0.082 J 
(4) Empirical Eq. (16 0.214 0.151 0.107 0.087 0.075 0.067 
b- | 


The experimental data on the lift of flat planing surfaces (ol 
tained by Sottorf, Sambraus, Shoemaker, and Locke) have been 
summarized in the empirical formula by Korvin-Kroukovsky, 
Savitsky, and Lehman in reference 3, 


Cr = r![0.0120\'/2 + 0.0095 (A/Cr)? l4 
where 
Cr = lift/(1/2)pb? V2 
b = beam 
A = ratioof wetted length to beam = //b 
Cy = speed coefficient, defined as V/V gb 
T = angle of trim with respect to the undisturbed water 


surface, in degrees 


Eq. (13) applies to the weightless fluid. In the experimental 
Eq. (14), this condition is approached at high values of Cy, at 
which the second term can be neglected. Physically, this corre- 
sponds to the range of speeds at which gravity forces can be neg- 
lected in comparison with dynamic forces. In Eq. (14), Ci 
(1/2)pb?V?, which should be changed for compat 
(13) to C. = L/(1/2)pbl1V2. With these changes, 


Eq. (14) becomes 


defined as L 
son with Eq 
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and 


dCi/dr = 0.013217™'/d 7? (16 


The data obtained from these empirical relations representing 
the experimental data are listed on lines 2 and 4 of Table 2 for 
comparison with the analytical data obtained by Eqs. (12) and 
(13). While the agreement is very good, it appears that the 
analytical formulas underestimate the last (nonlinear) term in 
Eq. (13) and thereby underestimate the lift at the high length 


beam ratios 
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Potential Flow Around a Rotating, Advancing 
Cylindrical Blade 


L. E. Fogarty and W. R. Sears 

Graduate Schoo! of Aeronautical Engineering, Cornell University 
Ithaca, N.Y 

June 14, 1950 


I A RECENT CONTRIBUTION,! Sears has treated the problem of 
an infinitely long cylinder rotating at constant angular velocity 
w about an axis normal to it in a perfect incompressible fluid 
otherwise at rest 

Using axes rotating with the blade, Sears measures =z along the 
axis of rotation and y outward along the blade from the center 
of rotation and takes x to form the third member of a right- 
handed Cartesian system. For y > 0 the cylinder moves in the 
—x direction. If u’, 7’, w’ are the velocity components meas 


ured with respect to the rotating system, then Sears’ results are 


u’ = wy¢iz(X, 2) / 
v’ = wlyi(x, 2) — 2x] > 1) 
w’ = wy¢iz(x, 2) \ 


=) denotes the potential for plane steady flow past the 


»< 


cylinder in a parallel stream of unit speed in the positive x direc 


where ¢; (x 


tion. 

It may be easily verified that these results can be extended to 
the case of a cylindrical blade advancing in the =z direction as it 
rotates about the z axis. Suppose the velocity of advance is WW’; 
then, keeping the coordinate system fixed in the blade, the ex 


tended results are 


u’ = WVPiz(X, 2) + W oz ( _# s) / 
= wly,(x, 2) — 2x] > (2) 
w’ = wygi2(x, 5) + Wee2(x, 2) \ 


where ¢o(x, 2) is the potential for plane steady flow past the cylin 
der in a parallel stream of unit speed in the negative =z direction 


<y 
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Thus, for a cylindrical blade advancing like a propeller, the 
flow components uw’ and w’ are the same as for two-dimensional 
motiou at the local relative speed and incidence. The spanwise 
flow, however, depends only on the potential distribution of a 
part of this two-dimensional flow. 

These results need not be restricted to flow without circulation 
as Sears implied in his note. For example, if the blade profile 
is that of a cambered airfoil, ¢)(x, 2) and go(x, =) each include a 
definite amount of circulation, determined by the Kutta-Jou- 
kowski condition. Thus, the total circulation about the rotating 
blade, to satisfy the trailing-edge condition in the superimposed 
flow, according to Eqs. (2), must be of the form 


r(y) = wy, + WP. 


where I’; and Il, are constants determined by the profile 

The presence of this variable circulation along the blade does 
not invalidate any step of Sears’ derivation,' but it implies that 
there is a helical trailing-vortex sheet of uniform strength behind 
the blade. We have not accounted for induced effects of such 
a sheet in satisfying the boundary condition at the blade. If the 
blade is symmetrical, T; is zero, the circulation is constant along 
the blade, there is no trailing wake, and all the conditions of the 


problem are satisfied 


REFERENCE 
Sears, W. R., Potential Flow Around a Rotating Cylindrical Blade, Read 
ers’ Forum, Journal of the Aeronautical Sciences, Vol. 17, No. 3, p. 183 
March, 1950 


Application of Biezeno-Koch Method to 
Compressible Fluid Flow Problems 


Chi-Teh Wang and Pei-Chi Chou 
Associate Professor and Graduate Student, respectively, 
Daniel Guggenheim School of Aeronautics, New York University 


June 16, 1950 


IT PREVIOUS PAPERS,' * the senior author and his associates 
have applied the variational methods to compressible fluid 
flow problems following the Rayleigh-Ritz and Galerkin pro 
cedures. For all the numerical examples worked out, very good 
agreement was obtained between the results computed by the 
variational methods and those by other approximate methods 
In this note, a method that was originally introduced by Biezeno 
and Koch® for the solution of plate and beam problems is now 
shown to give good approximate solution also in compressible 
flow problems. As in the case of the Galerkin method, the 
Biezeno-Koch method does not require the formulation of any 
variational principle. Further, it has the advantage over the 
Galerkin method in that the numerical work is simpler to carry 


out. 


BIEZENO-KocH METHOD 


The Biezeno-Koch method is to be carried out as follows 


First, assume a series of functions, containing undetermined 
parameters, which satisfy the boundary conditions. Next, sub 
stitute the assumed series into the governing differential equa 
tion. If the assumed series is the exact solution, the differential 
equation becomes identically zero upon substitution. But -if 
the assumed series is not the exact solution, as is usually the 
case, the resulting function is not zero and may be interpreted 


as an error function. In the so-called ‘‘collocation method,’’ one 


tries to minimize the error function by requiring it to vanish at 
N points of the region, where N is the number of parameters 
used. For two-dimensional problems, especially in the case of 
flow past an arbitrary body, where the domain is infinite, it 
appears that unless many points are taken the approximation 
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Mach 
Number 


Ay ay 
Rayleigh- 
Ritz* 
0.001031 
0.008583 
0.08114 
0.08307 

0.1979 


Biezeno- 
Koch 
0.001840 
0.015778 
0.060861 
0.18235 
+ 


Galerkin Koch 
0.001032 2.0123 
0.00859 2.0526 
0.03127 2.1352 
0.08360 2.3039 
0.19878 , 


*. 9 


+t No real solution exists (see reference 4). 


cannot be good. If so, the resulting simultaneous equations 
will be so many that it becomes impractical to solve them. In 
the Biezeno-Koch method, if N parameters are taken, the domain 
R is divided into V regions R;, and one minimizes the error func- 
tion by requiring the mean error over each of the N regions to 
be zero, thus giving .V equations for the determination of the N 
parameters. Since the mean errors in many regions are required 
to be zero rather than the errors at many specific points, it may 
be expected that the Biezeno-Koch method will give better 
| the 
numerical work. 

It was pointed out by Courant in the discussion of reference 10 


results than collocation method for the same amount of 


that the method of Biezeno and Koch is actually a special case 
Let 
The problem of 


of the variational method. This may be seen as follows. 
the governing differential equation be L = 0. 
taking the first variation of the variational integral J and the 


condition 6J = 0 leads to the equation 


if. Ln(x,y) dx dy = 0 


where 6 denotes the first variation, 7 is the arbitrary variation, 
and x,y are Cartesian coordinates. If one expands the arbitrary 
variation 7(x,y) in terms of the functions n;(x,y), which are de- 
fined by 

\ Lin R; 


i(x,y) = , 
” - { O elsewhere in R 


then one obtains exactly the Biezeno-Koch method. 


NUMERICAL EXAMPLE 


flow past a 
circular cylinder without circulation has been carried out. The 


fundamental differential equation of a compressible fluid in polar 


As an example, the problem of compressible 


coordinates is as follows: 


_ *) Pay 4 


r? r? 


6" or 9 dra 
, a <orh4 e 
r ial 


2 


(1) 


where ¢ is the velocity potential, aj is the velocity of sound at 
infinity, U is the undisturbed velocity, and the subscripts in- 


dicate partial differentiation. The boundary conditions are 
(¢,? + 1 52)' °= Uatr=o | 
@ = Oatr=1 \ 
The radius of the cylinder is taken to be unity. ¢@ may be 
assumed to be as follows: 


u(r+ ') cos 0 + =, = Amn X 
r 


m=zin=] 


1 l 
mm “4 9)ym + | cos n8 (3) 
nr (m + Z)r 


og = 
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, —e 
Gmez./ U- —__] 
, ; Rayleigh. 
Galerkin Jansen 
One Four (third ap. 
parameter parameters PrOximatiog) 
2.00688 > 2.01193 
2.0286 2.0524 2.02867 ; 2.05128 
2.0692 2.1364 2.06949 2.1314] 
2.1385 2.3336 2.13933 20 de 2. 28361 
2.2639 a 2.26504 2.57071 


5, 


Rayleigh-Ritz* 
One Six 
parameter parameters 

2.0069 2.0120 


Note that expression (3) satisfies the boundary conditions (¥ 
Substituting Eq. (3) (1) and calling the 
resulting error function 


exactly. into Eq. 
en = Ligon) 


where the subscript .V indicates the number of parameters taken, 
the Biezeno-Koch equations for determining the parameteg 


Sf. L(én)r dr dé = 0, +=1,2...N (4) 


Since the flow in this case is symmetrical with respect to both 
horizontal and vertical axes, only one-quarter of the domaig 
Taking one parameter A, only, Eq 


A mn are 


needs to be considered. 


(ey x eA =— 
J 1 J, = 0 


Carrying out the integration, we obtain 


(4) becomes 


~ L(g,)r dr do = 0 (5) 


0.0184404 1,3 + 0.389775 UA)? + 
(1.909905 U2 — 0.888889a0?) A + 1.6U% = 0 @& 


where y is taken to be 1.405. 

Solving at various Mach Numbers of the undisturbed flow 
My, Au/ao and the maximum velocity at the surface of the bod 
The agreement between the result§ 
methods 


are tabulated in Table 1. 
the 
appears to be very good. 


obtained by Biezeno-Koch method and other 
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